text:
#!/usr/bin/env julia
using Distributed
using MIToS.Utils.Scripts
using MIToS

Args = parse_commandline(
    ["--format", "-f"],
    Dict(
        :help => "Format of the MSA: Stockholm, Raw or FASTA",
        :arg_type => String,
        :default => "Stockholm"
    ),
    ["--exact", "-e"],
    Dict(
        :help => "If it's true, the mean percent identity is exact (using all the pairwise comparisons).",
        :action => :store_true
    ),
    description="""
    Creates an *.description.csv from a Stockholm file with: the number of columns, sequences, clusters after Hobohm clustering at 62% identity and mean percent identity.
    Also the mean, standard deviation and quantiles of: sequence coverage of the MSA, gap percentage.
    """,
    output=".description.csv",
	mitos_version=loadedversion(MIToS)
    # ----------------------------------------------------------------------------
    )

set_parallel(Args["parallel"])

@everywhere begin

    const args = remotecall_fetch(()->Args,1)

    import MIToS.Utils.Scripts: script

    # TO DO ----------------------------------------------------------------------
    using Pkg
    using Dates
    using Statistics
    using MIToS.Utils.Scripts
    using MIToS.MSA

    function _describe(fh_out, input, format)

    end
    # ----------------------------------------------------------------------------

    function script(input::Union{Base.LibuvStream,  AbstractString},
                    args,
                    fh_out::Union{Base.LibuvStream, IO})
        # TO DO ------------------------------------------------------------------
		println(fh_out, "# MIToS ", loadedversion(MIToS), " MSADescription.jl ", now())
        println(fh_out, "# used arguments:")

        for (key, value) in args
            println(fh_out, "# \t", key, "\t\t", value)
        end

        form = string(args["format"])

        if form == "Stockholm"
            aln = readorparse(input, Stockholm)
        elseif form == "FASTA"
            aln = readorparse(input, FASTA)
        elseif form == "Raw"
            aln = readorparse(input, Raw)
        else
            throw(ErrorException("--format should be Stockholm, Raw or FASTA."))
        end

        println(fh_out, input, ",", "sequences", ",", "number", ",", "", ",", size(aln, 1))
        println(fh_out, input, ",", "columns",   ",", "number", ",", "", ",", size(aln, 2))

        cov = vec(coverage(aln));
        qcv = quantile(cov, [0., .25, .5, .75, 1.])

        println(fh_out, input, ",", "coverage", ",", "quantile", ",", "0.00", ",", qcv[1])
        println(fh_out, input, ",", "coverage", ",", "quantile", ",", "0.25", ",", qcv[2])
        println(fh_out, input, ",", "coverage", ",", "quantile", ",", "0.50", ",", qcv[3])
        println(fh_out, input, ",", "coverage", ",", "quantile", ",", "0.75", ",", qcv[4])

        println(fh_out, input, ",", "coverage",   ",", "mean", ",", "", ",", mean(cov))
        println(fh_out, input, ",", "coverage",   ",", "std",  ",", "", ",", std(cov))

        println(fh_out, input, ",", "percentidentity",   ",", "mean", ",", "", ",", meanpercentidentity(aln, exact=args["exact"]))

        gap = vec(gapfraction(aln, 1));
        qgp = quantile(gap, [0., .25, .5, .75, 1.])

        println(fh_out, input, ",", "gapfraction", ",", "quantile", ",", "0.00", ",", qgp[1])
        println(fh_out, input, ",", "gapfraction", ",", "quantile", ",", "0.25", ",", qgp[2])
        println(fh_out, input, ",", "gapfraction", ",", "quantile", ",", "0.50", ",", qgp[3])
        println(fh_out, input, ",", "gapfraction", ",", "quantile", ",", "0.75", ",", qgp[4])

        println(fh_out, input, ",", "gapfraction",   ",", "mean", ",", "", ",", mean(gap))
        println(fh_out, input, ",", "gapfraction",   ",", "std",  ",", "", ",", std(gap))

        hob = hobohmI(aln, 62);
        ncu = nclusters(hob)

        println(fh_out, input, ",", "clusters",   ",", "number", ",", "", ",", ncu)
        # ------------------------------------------------------------------------
    end

end

runscript(args)


text:
mutable struct simulationTimer
  initialTime::Float64
  wallTime::Float64
  timeStepStartTime::Float64
  simulationTime::Float64
  finalTime::Float64
  stepIndex::Int64
  printIndex::Int64
  printSkip::Int64
  Δt::Float64
end

function exportVTK!(timer::simulationTimer,p::precipParam,soln::precipSolution,mesh::precipMesh)
  # scalar data
  sData = ScalarData(soln.ψA.u,soln.ψB.u,soln.ψC.u,soln.θs.u,soln.θf.u)
  sNames = ScalarNames(p.ψAname,p.ψBname,p.ψCname,p.θsname,p.θfname)
  path = Path(p.SaveFolder,string(p.SaveName,"_scalars_",timer.printIndex))
  vtksave(mesh.sMesh,sData,sNames,path;time=timer.simulationTime)

  # fluids data
  sData = ScalarData(soln.fluid.p)
  sNames = ScalarNames("pressure")
  vData = VectorData([soln.fluid.u,soln.fluid.v])
  vNames = VectorNames("darcy_velocity")
  path = Path(p.SaveFolder,string(p.SaveName,"_fluid_",timer.printIndex))
  vtksave(mesh.fMesh,sData,sNames,vData,vNames,path;time=timer.simulationTime)
end

function doPrinting!(timer::simulationTimer,p::precipParam,mesh::precipMesh,soln::precipSolution)
  if p.printBool
    printProgressBar(timer)
    exportVTK!(timer,p,soln,mesh)
  end
  updatePrintIndex!(timer)
end

function initializeTimer(p::precipParam)
    return simulationTimer(time(),0.0,0.0,eps(),p.Tfinal,1,1,p.printSkipIndex,p.Δt)
end

function updateWallTime!(timer::simulationTimer)
    timer.wallTime = time() - timer.initialTime
end
function updateSimulationTime!(timer::simulationTimer)
    timer.simulationTime += timer.Δt
end
function updateStepIndex!(timer::simulationTimer)
    timer.stepIndex += 1
end
function updatePrintIndex!(timer::simulationTimer)
    timer.printIndex += 1
end
function intializeTimeSteppingTimer!(timer::simulationTimer)
  timer.timeStepStartTime::Float64
end

function isPrintIndex(timer::simulationTimer)
    if mod(timer.stepIndex,timer.printSkip)==0
        return true
    end
    return false
end

function printSimulationFinishedMessage(timer::simulationTimer)
  println()
  println("  Simulation done ... took ",sToHMS(timer.wallTime))
end

function printProgressBar(timer::simulationTimer)
  percentage = Int(ceil(timer.simulationTime/timer.finalTime*100))
  progressBarDisplay(percentage,timer)
end

function progressBarDisplay(percentage,timer::simulationTimer)
  TERMINALLENGTH = 80

  percent = string("\r  [",percentage,"%]")
  skip1  = 11-length(percent)
  exportFigure = string("exported fig. ",timer.printIndex)
  skip2  = 31-(length(percent)+skip1+length(exportFigure))
  updateWallTime!(timer)
  timeTaken = string("wall time: ",sToHMS(timer.wallTime),)
  skip3  = 48-(length(percent)+skip1+length(exportFigure)+skip2+length(timeTaken))
  tremain = (timer.finalTime - timer.simulationTime)*(timer.wallTime - timer.timeStepStartTime)/timer.simulationTime
  timeRemaining1 = string("remaining: ")
  timeRemaining2 = string(sToHMS(tremain),"\r")

  iterChar(TERMINALLENGTH," ")

  print(percent)

  space()
  iterChar(skip1,"-")
  space()
  print(exportFigure)

  space()
  iterChar(skip2,"-")
  space()
  print(timeTaken)

  space()
  iterChar(skip3,"-")
  space()
  print(timeRemaining1)
  if timer.printIndex > 4
    print(timeRemaining2)
  else
    print("...\r")
  end
end

function iterChar(iterations::Int,character::String)
  for i=1:iterations
    print(character)
  end
end

function space()
    print(" ")
end

function doInitialPrinting(p::precipParam)
  #if p.printBool
    println("  running ",returnRunIDName(p))
    println("  data being saved in: ",p.SaveFolder)

    run(`mkdir -p $(p.SaveFolder)`)
    exportParameters(p)
  #end
end

function returnRunIDName(p::precipParam)
  if p.runID == 1
    return "main program"
  elseif p.runID == 2
    return "pure diffusion one chemical test"
  elseif p.runID == 3
    return "pure advection one chemical test"
  elseif p.runID == 4
    return "diffusion + advection one chemical test"
  elseif p.runID == 5
    return "multiphase brinkman no chemical test"
  elseif p.runID == 6
    return "diffusion + reaction two chemical test"
  elseif p.runID == 7
    return "pure reaction test"
  else
    error("runID not recognized")
  end
end

"""
  sToHMS(s::Float64)

Converts Float64 in seconds to HH:MM:SS string format
"""
function sToHMS(s)
  minutes = 0.0
  hours   = 0.0
  timestr = ""

  # compute times
  hours   = Int(floor(s/3600))
  minutes = Int(floor(s/60 - 60*hours))
  seconds = Int(round(s - 3600*hours - 60*minutes,digits=0))

  hours < 10   ? hstr = "0" : hstr = ""
  minutes < 10 ? mstr = "0" : mstr = ""
  seconds < 10 ? sstr = "0" : sstr = ""

  if hours > 0
    timestr=string(timestr,hstr,hours,"h:")
  end
  if (minutes>0) || (hours > 0)
    timestr=string(timestr,mstr,minutes,"m:")
  end

  timestr = string(timestr,sstr,seconds,"s")

  return timestr
end

function exportParameters(p::precipParam)
  dict = type2dict(p)

  filename = string(p.SaveFolder,"/",p.SaveName,"_parameters.txt")
  open(filename, "w") do io
    for (key,value) in dict
      str=string(key,": ",value,"\n")
      write(io,str)
    end
  end

  nothing
end

ext: jl


text:
push!(LOAD_PATH,"../src/")

using Documenter
using Scruff

makedocs(
    sitename = "Scruff.jl",
    modules = [Scruff, Scruff.SFuncs, Scruff.Algorithms, Scruff.Models,
        Scruff.Operators, Scruff.Utils, Scruff.RTUtils],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Avi Pfeffer, Michael Harridon, Joseph Campolongo, Sanja Cvijic, and contributors.",
    pages = [
        "Getting Started" => "index.md",
	    "Tutorial" => "tutorial/tutorial.md",
        "Examples" => "tutorial/examples.md",
        "Library" => Any[
            "Core" => "lib/core.md",
            "Stochastic Functions" => "lib/sfuncs.md",
            "Operators" => "lib/operators.md",
            "Models" => "lib/models.md",
            "Algorithms" => "lib/algorithms.md",
            "Utilities" => "lib/utilities.md",
            "Runtime Utilities" => "lib/rtutils.md"
        ]
    ]
)

deploydocs(
    repo="github.com/p2t2/Scruff.jl",
    devbranch = "develop")


text:
module Grids

export RegularCartesianGrid

using Oceananigans

include("grid_utils.jl")
include("regular_cartesian_grid.jl")

end


text:
export ZDT, ZDT1, ZDT2, ZDT3, ZDT4, ZDT6;

zdt_ref_string = "Zitzler et al. - 2000 - Comparison of Multiobjective Evolutionary Algorithms";
@doc """
    abstract type ZDT <: MOP end

Abstract super type for the class of ZDT problems as defined in
$(zdt_ref_string)

These are biobjective problems with a variable number of decision variables
`n_vars >= 2`.
They share a common structure
```math
\\begin{aligned}
    f_1(\\mathbf{x}) &= f( x_1 ),\\\\
    f_2(\\mathbf{x}) &= g( x_2, …, x_n ) · h( f_1( x_1 ), g( x_2, …, x_n ) )
\\end{aligned}
```
"""
abstract type ZDT <: MOP end;

# some ipmlementations that all problems share
num_vars( mop :: Z where Z<:ZDT ) = mop.n_vars; # make sure below that every implementation has field n_vars
num_objectives( mop :: Z where Z<:ZDT ) = 2;

get_pareto_set( mop :: Z where Z<:ZDT ) = SamplingFunction( mop, :ParetoSet, [:regular, :random] );
get_pareto_front( mop :: Z where Z<:ZDT ) = SamplingFunction( mop, :ParetoFront, [:regular, :random] );

function get_objectives( mop :: Z where Z<:ZDT )
    f1, g, h = zdt_helper_functions( mop )
    f2(x :: Vector{R} where R<:Real) = g(x) * h(f1(x), g(x))
    return [f1; f2]
end

########## ZDT 1 ##############
@doc """
    ZDT1( n_vars :: Int )

Bi-Objective problem ZDT1 with convex Pareto Front as given in
$(zdt_ref_string).

The objectives are given as
```math
\\begin{aligned}
    f_1( x_1 ) &= x_1,\\\\
    f_2( x_1, …, x_n) &= g( x_2, …, x_n ) · h( f_1( x_1 ), g( x_2, …, x_n ) ), \\\\
    g( x_2, …, x_n ) &= 1 + (9/(n-1)) · Σ_i x_i, \\\\
    h( f1, g ) &= 1 - \\sqrt{f_1/g}.
\\end{aligned}
```
"""
@with_kw struct ZDT1 <: ZDT
    n_vars :: Int = 2;
end

constraints( mop :: ZDT1 ) = Box( zeros(mop.n_vars), ones( mop.n_vars) );

function zdt_helper_functions( mop :: ZDT1 )
    h(f1,g) = 1 - sqrt( f1 / g );
    g(x) = 1 + 9 / (num_vars(mop) - 1) * sum(x[2:end]);
    f1(x :: Vector{R} where R<:Real) = x[1]
    return [f1, g, h]
end

function get_ideal_point(mop::ZDT1)
    # min of f1 = x is mop.lb[1] = 0
    # min of g is 1, min of h is 1 - sqrt(ub[1]) = 0
    return zeros(2);
end

########## ZDT 2 ##############
@doc """
    ZDT2( n_vars :: Int )

Bi-Objective problem ZDT2 with **non**convex Pareto Front as given in
$(zdt_ref_string).

The objectives are given as
```math
\\begin{aligned}
    f_1( x_1 ) &= x_1,\\\\
    f_2( x_1, …, x_n) &= g( x_2, …, x_n ) · h( f_1( x_1 ), g( x_2, …, x_n ) ), \\\\
    g( x_2, …, x_n ) &= 1 + (9/(n-1)) · Σ_i x_i, \\\\
    h( f1, g ) &= 1 - (f_1/g)².
\\end{aligned}
```
"""
@with_kw struct ZDT2 <: ZDT
    n_vars :: Int = 2;
end

constraints( mop :: ZDT2 ) = Box( zeros(mop.n_vars), ones( mop.n_vars) );

function zdt_helper_functions(mop::ZDT2)
    h(f1,g) = 1 - ( f1 / g )^2;
    g(x) = 1 + 9 / (num_vars(mop) - 1) * sum(x[2:end]);
    f1(x :: Vector{R} where R<:Real) = x[1]
    return [f1, g, h]
end

get_ideal_point(mop::ZDT2) = zeros(2);

########## ZDT 3 ##############
@doc """
    ZDT3( n_vars :: Int )

Bi-Objective problem ZDT3 with discrete Pareto Front as given in
$(zdt_ref_string).

The objectives are given as
```math
\\begin{aligned}
    f_1( x_1 ) &= x_1,\\\\
    f_2( x_1, …, x_n) &= g( x_2, …, x_n ) · h( f_1( x_1 ), g( x_2, …, x_n ) ), \\\\
    g( x_2, …, x_n ) &= 1 + (9/(n-1)) · Σ_i x_i, \\\\
    h( f1, g ) &= 1 - \\sqrt{f_1/g}- (f_1/g)·\\sin( 10 π f_1 )
\\end{aligned}
```
"""
@with_kw struct ZDT3 <: ZDT
    n_vars :: Int = 2;
end

constraints( mop :: ZDT3 ) = Box( zeros(mop.n_vars), ones( mop.n_vars) );

function zdt_helper_functions(mop::ZDT3)
    h(f1,g) = 1 - sqrt( f1 / g ) - (f1 / g) * sin(10 * pi * f1);
    g(x) = 1 + 9 / (num_vars(mop) - 1) * sum(x[2:end]);
    f1(x :: Vector{R} where R<:Real) = x[1]

    return [f1, g, h]
end
get_ideal_point(mop::ZDT3) = zeros(2);

########## ZDT 4 ##############
@doc """
    ZDT4( n_vars :: Int )

Bi-Objective problem ZDT4 with many local Pareto Fronts as given in
$(zdt_ref_string).

The multimodality is due to ``g`` being a version of Rastrigin's function.

The objectives are given as
```math
\\begin{aligned}
    f_1( x_1 ) &= x_1,\\\\
    f_2( x_1, …, x_n) &= g( x_2, …, x_n ) · h( f_1( x_1 ), g( x_2, …, x_n ) ), \\\\
    g( x_2, …, x_n ) &= 1 + (9/(n-1)) · Σ_i ( x_i^2 - \\cos( 4 π x_i ) ), \\\\
    h( f1, g ) &= 1 - \\sqrt{f_1/g}
\\end{aligned}
```
"""
@with_kw struct ZDT4 <: ZDT
    n_vars :: Int = 2;
end

constraints( mop :: ZDT4 ) = Box(
    [ 0; fill(-5.0, mop.n_vars - 1) ],
    [ 1.0; fill(5.0, mop.n_vars - 1) ]
);

function zdt_helper_functions(mop::ZDT4)
    h(f1,g) = 1 - sqrt( f1 / g )
    g(x) = 1 + 10 * ( num_vars(mop) -1 ) + sum( x[2:end].^2 -10*cos.( 4*pi .* x[2:end]) )
    f1(x :: Vector{R} where R<:Real) = x[1]

    return [f1, g, h]
end
get_ideal_point(mop::ZDT4) = zeros(2);

get_critical_set( mop :: Z where Z<:ZDT ) = SamplingFunction( mop, :ParetoCriticalSet, [:regular, :random] );

########## ZDT 5 ##############
# not implemented yet due to it having binary variables

########## ZDT 6 ##############
@doc """
    ZDT6( n_vars :: Int )

Bi-Objective problem ZDT6 with biased Pareto Front as given in
$(zdt_ref_string).

The objectives are given as
```math
\\begin{aligned}
    f_1( x_1 ) &= 1 - \\exp(-4 x_1 ) \\sin^6( 6 π x_1) ,\\\\
    f_2( x_1, …, x_n) &= g( x_2, …, x_n ) · h( f_1( x_1 ), g( x_2, …, x_n ) ), \\\\
    g( x_2, …, x_n ) &= 1 + 9 · Σ_i ( x_i / (n-1) )^{1/4}, \\\\
    h( f1, g ) &= 1 - (f_1/g)²
\\end{aligned}
```
"""
@with_kw struct ZDT6 <: ZDT
    n_vars :: Int = 2;
end

constraints( mop :: ZDT6 ) = Box( zeros(mop.n_vars), ones( mop.n_vars) );

function zdt_helper_functions(mop::ZDT6)
    h(f1,g) = 1 - ( f1 / g )^2;
    g(x) = 1 + 9 * ( sum(x[2:end])/(num_vars(mop)-1) )^0.25;
    f1(x :: Vector{R} where R<:Real) = x[1] - exp( -4*x[1] ) * (sin(6*pi*x[1]))^6;

    return [f1, g, h]
end
get_ideal_point( mop::ZDT6 ) = zeros(2);

####### Point Sampling Functions #########################

#helper
function zdt_x_range(mop :: Z where Z<:ZDT, method :: Symbol, num_points::Int,
    x1 :: Real, x2 :: Union{Nothing,Vector{R} where R<:Real} = nothing )

    if isnothing(x2)
        x2 = zeros( num_vars(mop) -1 );
    end

    @assert length(x2) == num_vars(mop) - 1;

    if method == :regular
        x_range = num_points == 1 ? [x1,] : range(x1, 1.0; length = num_points );
    elseif method == :random
        x_range = [ x1 + (1-x1) * rand() for i = 1 : num_points ]
    end
    return [ [x; x2] for x in x_range ]
end

@doc """
    get_points( sample_func, zdt, Val(:ParetoSet), N; method = :regular )

Return a list of points ``x ∈ ℝⁿ`` where ``x₁ = … = xₙ = 0``.
Methods are `:regular` or `:random`.
"""
function get_points( sample_func :: SamplingFunction, mop :: Union{ZDT1,ZDT2,ZDT4},
        ::Val{:ParetoSet}, num_points :: Int; method = :regular )
    zdt_x_range(mop, method, num_points, 0.0)
end

function get_points( sample_func :: SamplingFunction, mop :: Union{ZDT6},
        ::Val{:ParetoSet}, num_points :: Int; method = :regular )
    # see
    # https://sop.tik.ee.ethz.ch/download/supplementary/testproblems/zdt6/index.php
    x1 = .2807753191;
    zdt_x_range( mop, method, num_points, x1)
end

# ZDT 3 is a special snowflake …
function sample_from_range( arr, N, ::Val{:regular} )
    [ arr[1] + (arr[2] - arr[1]) * x for x in range(0, 1;length = N) ]
end

function sample_from_range( arr, N, ::Val{:random} )
    [ arr[1] + (arr[2] - arr[1]) * rand() for i = 1 : N ]
end

function get_points( :: SamplingFunction, mop :: ZDT3, ::Val{:ParetoSet},
        num_points :: Int; method = :regular )
    # see
    # https://sop.tik.ee.ethz.ch/download/supplementary/testproblems/zdt3/index.php

    chunk_size = num_points ÷ 5;

    x_ranges = [
        [0, 0.0830015349],
            [[1e-10, 0.0] .+ v
            for
            v ∈ [
                [0.1822287280, 0.2577623634],
                [0.4093136748, 0.4538821041],
                [0.6183967944, 0.6525117038],
                [0.8233317983, 0.8518328654],
            ]]...
    ];

    x_range = Float64[];
    for i = 1 : 5
        if i == 1
            N = chunk_size + (num_points % chunk_size)
        else
            N = chunk_size
        end
        push!( x_range, sample_from_range( x_ranges[i], N, Val(method) )... )
    end
    return [ [x; zeros( num_vars(mop) -1 )] for x in x_range ]
end

function get_points( :: SamplingFunction, mop :: ZDT4, ::Val{:ParetoCriticalSet},
        num_points :: Int; x0 = nothing, method = :regular )
    #@eval using NLopt;

    if isnothing(x0)
        x̂ = zeros(mop.n_vars - 1) # global optimum;
    else
        if length(x0) == mop.n_vars
            x0 = x0[2:end];
        end

        if mop.n_vars == 2 && isa( x0, Real )
            x0 = [x0,]
        end

        _, g, _ = zdt_helper_functions( mop );
        # gradient via autodiff
        # g does not depend on x1 = 0.0
        ∇g = x -> ForwardDiff.gradient( g, [0.0; x] )[2:end];

        opt_func = function( x::Vector, grad::Vector)
            if length(grad) > 0
                grad[:] = ∇g(x)
            end
            return g(x)
        end

        box = constraints( mop );

        opt = NLopt.Opt( :LD_MMA, length(x0) );
        opt.lower_bounds = box.lb[2:end];
        opt.upper_bounds = box.ub[2:end];

        opt.xtol_rel = 1e-5;
        opt.maxeval = 200 * length(x0);

        println("Starting search for local minimum of g for ZDT4 with x0:")

        opt.min_objective = opt_func;
        _, x̂, ret =  NLopt.optimize( opt, x0 );
        #@show x̂
        println("Finished with $(opt.numevals) eval(s).")
    end

    return zdt_x_range( mop, method, num_points, 0.0, x̂ );
end

get_critical_front( mop :: Z where Z<:ZDT ) = SamplingFunction( mop, :ParetoCriticalFront, [:regular, :random] );


text:
"""
Creates a python file that can be run in Blender to generate a
representation of an object created by the CGinJulia package

'''Bl_Export(filename::String, points::AbstractVector{Point3H}, faces::AbstractVector{Real})'''
"""
function Bl_Export(filename::String,points,faces)

io = open(filename, "w")

Code::String = """bl_info = {
    "name": "Convex Hull Generation in Julia",
    "author": "Jackson Meade (NCSSM) and Jack Snoeyink (UNC Chapel Hill)",
    "version": (1, 0),
    "blender": (2, 80, 0),
    "location": "View3D > Add > Mesh > CGinJulia Convex Hull",
    "description": "Adds a new Mesh Object",
    "warning": "Please make sure you are familiar with the CGinJulia package before using this EXPERIMENTAL add-on",
    "wiki_url": "github.com/snoeyink/CGinJulia",
    "category": "Add Convex Hull",
}

import bpy
from bpy.types import Operator
from bpy.props import FloatVectorProperty
from bpy_extras.object_utils import AddObjectHelper, object_data_add
from mathutils import Vector

def add_object(self, context):
    scale_x = self.scale.x
    scale_y = self.scale.y
    scale_z = self.scale.z

    verts = ["""

for i = 1:length(points)
    Code = string(Code,string("""Vector((""",points[i][2],""" * scale_x""",""",""",points[i][3],""" * scale_y""",""",""",points[i][4],""" * scale_z""",""")),""","""
    """))
end

Code = string(Code,"""]
""","""    edges = []
""","""    faces = [""")

for a in faces
    Code = string(Code,string("""[""",string(a[1]),""",""",string(a[2]),""",""",string(a[3]),"""],"""))
end

Code = string(Code,"""]
""","""    mesh = bpy.data.meshes.new(name="Convex Hull of """,filename,"""")
    mesh.from_pydata(verts, edges, faces)
    # useful for development when the mesh may be invalid.
    # mesh.validate(verbose=True)
    object_data_add(context, mesh, operator=self)

class OBJECT_OT_add_object(Operator, AddObjectHelper):
    """,'\u0022','\u0022','\u0022',"""Generate the Mesh Object from Points and CGinJulia-Constructed Faces""",'\u0022','\u0022','\u0022',"""

    bl_idname = "mesh.add_object"
    bl_label = "Show Representation"
    bl_options = {'REGISTER', 'UNDO'}

    scale: FloatVectorProperty(
        name="scale",
        default=(1.0, 1.0, 1.0),
        subtype='TRANSLATION',
        description="scaling",
    )

    def execute(self, context):

        add_object(self, context)

        return {'FINISHED'}

# Registration

def add_object_button(self, context):
    self.layout.operator(
        OBJECT_OT_add_object.bl_idname,
        text="CGinJulia Convex Hull (""",filename,""")",
        icon='SURFACE_NSPHERE')

# This allows you to right click on a button and link to the manual
def add_object_manual_map():
    url_manual_prefix = "https://docs.blender.org/manual/en/dev/"
    url_manual_mapping = (
        ("bpy.ops.mesh.add_object", "editors/3dview/object"),
    )
    return url_manual_prefix, url_manual_mapping

def register():
    bpy.utils.register_class(OBJECT_OT_add_object)
    bpy.utils.register_manual_map(add_object_manual_map)
    bpy.types.VIEW3D_MT_mesh_add.append(add_object_button)

def unregister():
    bpy.utils.unregister_class(OBJECT_OT_add_object)
    bpy.utils.unregister_manual_map(add_object_manual_map)
    bpy.types.VIEW3D_MT_mesh_add.remove(add_object_button)

if __name__ == "__main__":
    register()
""")

write(io, Code)

close(io)

end

max_stars_repo_name: snoeyink/CGinJulia

max_issues_repo_name: snoeyink/CGinJulia

max_forks_repo_name: snoeyink/CGinJulia


text:
export unwrap_fun, wrap_fun, unwrap_head, wrap_head, unwrap_fcall, wrap_fcall
################################################################
"""
    head, body = unwrap_fun(fexpr)
    fcall, wherestack, body = unwrap_fun(fexpr,true)
    f, args, wherestack, body = unwrap_fun(fexpr, true, true)

Unwraps function expression.
"""
function unwrap_fun(expr::Expr)
    if expr.head in (:function, :(=))
        fexpr = expr
    elseif expr.head == :block
        fexpr = expr.args[2] # separate fexpr from block
    else
        error("Expression is not supported")
    end

    head = fexpr.args[1]
    body = fexpr.args[2]
    return head, body
end

function unwrap_fun(expr::Expr, should_unwrap_head::Bool)
    if expr.head in (:function, :(=))
        fexpr = expr
    elseif expr.head == :block
        fexpr = expr.args[2] # separate fexpr from block
    else
        error("Expression is not supported")
    end

    head = fexpr.args[1]
    fcall, wherestack = unwrap_head(head)
    body = fexpr.args[2]
    return fcall, wherestack, body
end

function unwrap_fun(expr::Expr, should_unwrap_head::Bool, should_unwrap_fcall::Bool)
    if expr.head in (:function, :(=))
        fexpr = expr
    elseif expr.head == :block
        fexpr = expr.args[2] # separate fexpr from block
    else
        error("Expression is not supported")
    end

    head = fexpr.args[1]
    fcall, wherestack = unwrap_head(head)
    f, args = unwrap_fcall(fcall)

    body = fexpr.args[2]
    return f, args, wherestack, body
end
################################################################
"""
    fexpr = wrap_fun(f, args, wherestack, body)
    fexpr = wrap_fun(fcall, wherestack, body)
    fexpr = wrap_fun(head, body)
    fexpr = wrap_fun(fexpr)

Returns a function definition expression
"""
function wrap_fun(f, args, wherestack, body)
    fcall = wrap_fcall(f, args)
    head =  wrap_head(fcall, wherestack)
    return Expr(:function, head, Expr(:block, body))
end

function wrap_fun(fcall, wherestack, body)
    head =  wrap_head(fcall, wherestack)
    return Expr(:function, head, Expr(:block, body))
end

function wrap_fun(head::Expr, body::Expr)
    return Expr(:function, head, Expr(:block, body))
end

function wrap_fun(fexpr::Expr)
    if fexpr.head in (:function, :(=))
        return fexpr
    elseif fexpr.head == :block
        fexpr = fexpr.args[2] # separate fexpr from block
        return fexpr
    else
        error("Expression is not supported")
    end
end

################################################################
function unwrap_head(head)
    wherestack = Any[]
    while head isa Expr && head.head == :where
        push!(wherestack, head.args[2])
        head = head.args[1]
    end
    fcall = head
    fcall, wherestack
end

function wrap_head(fcall, wherestack)
    for w in Iterators.reverse(wherestack)
        fcall = Expr(:where, fcall, w)
        # fcall = Expr(:where, fcall, esc(w))
    end
    head = fcall
    return head
end
################################################################
function unwrap_fcall(fcall::Expr)
    if !(fcall.head == :call)
        error("Expression is not supported")
    end
    f = fcall.args[1]
    args = fcall.args[2:end]
    return f, args
end

function wrap_fcall(f, args)
    fcall = :($f($((args)...)))
    return fcall
end
################################################################


text:
using DataDeps

register(DataDep("WikiGold",
    """
    Website: https://figshare.com/articles/Learning_multilingual_named_entity_recognition_from_Wikipedia/5462500

    WikiGold is a manually annotated corpus over a small sample of Wikipedia articles in CoNLL format (IOB).
    It contains 145 documents, 1696 sentences and 39152 tokens.
    The words have been labelled for each of the four CONLL-03 Named Enitty classes (LOC, MISC, ORG, PER).

    Please cite the following publication, if you are using the corpora:
    Dominic Balasuriya, Nicky Ringland, Joel Nothman, Tara Murphy, James R. Curran, 2009. Named Entity Recognition in Wikipedia, Proceedings of the 2009 Workshop on the People’s Web Meets NLP, ACL-IJCNLP 2009, pages 10–18.
    https://www.aclweb.org/anthology/papers/W/W09/W09-3302/
    """,
    "https://s3-eu-west-1.amazonaws.com/pfigshare-u-files/9446377/wikigold.conll.txt",
    "c797a64d0cf73ed058363f77671486bcfd413e70fda1726ad4a6ba624455225f";
    post_fetch_method = function(fn)
    end
))


text:
#' ---
#' title: Geometric Problems
#' ---

#' **Originally Contributed by**: Arpit Bhatia

#' These problems in this tutorial are drawn from Chapter 8 of the book
#' Convex Optimization by Boyd and Vandenberghe[[1]](#c1)

using JuMP
using Ipopt
# for plots
using Gadfly
using DataFrames

#' ## Euclidean Projection on a Hyperplane
#' For a given point $x_{0}$ and a set $C$, we refer to any point $z \in C$
#' which is closest to $x_{0}$ as a projection of $x_{0}$ on $C$.
#' The projection of a point $x_{0}$ on a hyperplane $C = \{x | a' \cdot x = b\}$ is given by

#' $$
#' \begin{align*}
#' \min && ||x - x_{0}|| \\
#' s.t. && a' \cdot x = b
#' \end{align*}
#' $$

x = rand(10)
a = rand(10)
b = rand()

projection = Model(with_optimizer(Ipopt.Optimizer, print_level=0))
@variable(projection, x0[1:10])
@objective(projection, Min, sum((x - x0) .* (x - x0))) # We minimize the square of the distance here
@constraint(projection, x0' * a == b)                  # Point must lie on the hyperplane

optimize!(projection)
@show objective_value(projection);
@show value.(x0);

#' ## Euclidean Distance Between Polyhedra
#' Given two polyhedra $C = \{x | A_{1} \cdot x \leq b1\}$ and $D = \{x | A_{2} \cdot x \leq b_{2}\}$,
#' the distance between them is the optimal value of the problem:

#' $$
#' \begin{align*}
#' \min && ||x - y|| \\
#' s.t. && A_{1} \cdot x \leq b_{1} \\
#' && A_{2} \cdot y \leq b_{2}
#' \end{align*}
#' $$

A_1 = rand(10, 10)
A_2 = rand(10, 10)
b_1 = rand(10)
b_2 = rand(10)

polyhedra_distance = Model(with_optimizer(Ipopt.Optimizer, print_level=0))
@variable(polyhedra_distance, x[1:10])                       # Point closest on the first polyhedron
@variable(polyhedra_distance, y[1:10])                       # Point closest on the second polyhedron
@objective(polyhedra_distance, Min, sum((x - y) .* (x - y))) # We minimize the square of the distance here as above
@constraint(polyhedra_distance, A_1 * x .<= b_1)             # Point x must lie on the first polyhedron
@constraint(polyhedra_distance, A_2 * y .<= b_2)             # Point y must lie on the second polyhedron

optimize!(polyhedra_distance)
@show objective_value(polyhedra_distance);

#' ## Linear Placement Problem
#' We have $N$ points in $\mathbb{R}^2$, and a list of pairs of points that must be connected by links.
#' The positions of some of the $N$ points are fixed; our task is to determine the positions of the remaining points,
#' i.e., to place the remaining points. The objective is to place the points so that the distance between the links is minimized,
#' i.e. our objective is:

#' $$
#' \begin{align*}
#' \min && \sum_{(i,j) \in A}||p_{i} - p_{j}||
#' \end{align*}
#' $$

fixed = [ 1   1  -1 -1    1   -1  -0.2  0.1;         # coordinates of fixed points
          1  -1  -1  1 -0.5 -0.2    -1    1]

M = size(fixed,2)                                    # number of fixed points
N = 6                                                # number of free points

A = [ 1  0  0 -1  0  0    0  0  0  0  0  0  0  0;    # Matrix on links
      1  0 -1  0  0  0    0  0  0  0  0  0  0  0;
      1  0  0  0 -1  0    0  0  0  0  0  0  0  0;
      1  0  0  0  0  0   -1  0  0  0  0  0  0  0;
      1  0  0  0  0  0    0 -1  0  0  0  0  0  0;
      1  0  0  0  0  0    0  0  0  0 -1  0  0  0;
      1  0  0  0  0  0    0  0  0  0  0  0  0 -1;
      0  1 -1  0  0  0    0  0  0  0  0  0  0  0;
      0  1  0 -1  0  0    0  0  0  0  0  0  0  0;
      0  1  0  0  0 -1    0  0  0  0  0  0  0  0;
      0  1  0  0  0  0    0 -1  0  0  0  0  0  0;
      0  1  0  0  0  0    0  0 -1  0  0  0  0  0;
      0  1  0  0  0  0    0  0  0  0  0  0 -1  0;
      0  0  1 -1  0  0    0  0  0  0  0  0  0  0;
      0  0  1  0  0  0    0 -1  0  0  0  0  0  0;
      0  0  1  0  0  0    0  0  0  0 -1  0  0  0;
      0  0  0  1 -1  0    0  0  0  0  0  0  0  0;
      0  0  0  1  0  0    0  0 -1  0  0  0  0  0;
      0  0  0  1  0  0    0  0  0 -1  0  0  0  0;
      0  0  0  1  0  0    0  0  0  0  0 -1  0  0;
      0  0  0  1  0  0    0  0  0  0  0 -1  0  0;
      0  0  0  0  1 -1    0  0  0  0  0  0  0  0;
      0  0  0  0  1  0   -1  0  0  0  0  0  0  0;
      0  0  0  0  1  0    0  0  0 -1  0  0  0  0;
      0  0  0  0  1  0    0  0  0  0  0  0  0 -1;
      0  0  0  0  0  1    0  0 -1  0  0  0  0  0;
      0  0  0  0  0  1    0  0  0  0 -1  0  0  0;]

placement = Model(with_optimizer(Ipopt.Optimizer, print_level=0))
@variable(placement, p[1:M + N, 1:2])                     # A variable array for the coordinates of each point
@constraint(placement, p[N + 1:N + M, :] .== fixed')      # We had a constraint for the fixed points
dist = A * p                                         # Matrix of differences between coordinates of 2 points with a link
@objective(placement, Min, sum(dist .* dist))            # We minimize the sum of the square of the distances

optimize!(placement)
@show value.(p);
@show objective_value(placement);

#+

# Plotting the points
df = DataFrame()
df.x = value.(p)[:,1]
df.y = value.(p)[:,2]
df.type = vcat(fill("Free points", N), fill("Fixed points", M))
plt = plot(df, x = "x", y = "y", color = "type", Geom.point)
draw(SVG(6inch, 6inch), plt)

#' ## Floor Planning
#' A floor planning problem consists of rectangles or boxes aligned with the axes which must be placed,
#' within some limits such that they do not overlap. The objective is usually to minimize the size
#' (e.g., area, volume, perimeter) of the bounding box, which is the smallest box that contains the boxes to be configured and placed.
#' We model this problem as follows:

#' We have N boxes $B_{1}, . . . , B_{N}$ that are to be configured and placed in a rectangle with width $W$ and height $H$,
#' and lower left corner at the position $(0, 0)$. The geometry and position of the $i$th box is specified by
#' its width $w_{i}$ and height $h_{i}$, and the coordinates $(x_{i}, y_{i})$ of its lower left corner.

#' The variables in the problem are $x_{i}, y_{i}, w_{i}, h_{i}$ for $i=1, \ldots,
#' N,$ and the width $W$ and height $H$ of the bounding rectangle. In all floor planning problems,
#' we require that the cells lie inside the bounding rectangle, $i . e .$

#' $$
#' x_{i} \geq 0, \quad y_{i} \geq 0, \quad x_{i}+w_{i} \leq W, \quad y_{i}+h_{i} \leq H, \quad i=1, \ldots, N
#' $$

#' We also require that the cells do not overlap, except possibly on their boundaries, i.e.

#' $$
#' x_{i}+w_{i} \leq x_{j} \quad \text{or} \quad x_{j}+w_{j} \leq x_{i} \quad \text{or} \quad y_{i}+h_{j} \leq y_{j} \quad \text{or} \quad y_{j}+h_{i} \leq y_{i}
#' $$

n = 5;

Amin = [                                        # We'll try this problem with 4 times with different minimum area constraints
100 100 100 100 100;
 20  50  80 150 200;
180  80  80  80  80;
 20 150  20 200 110]

r = 1

figs=[]

for i = 1:4
    A = Amin[i, :]

    floor_planning = Model(with_optimizer(Ipopt.Optimizer, print_level=0))

    @variables(floor_planning, begin
        x[1:n] >= r
        y[1:n] >= r
        w[1:n] >= 0
        h[1:n] >= 0
        W
        H
    end)

    @constraints(floor_planning, begin
        x[5] + w[5] + r <= W    # No rectangles at the right of Rectangle 5
        x[1] + w[1] + r <= x[3] # Rectangle 1 is at the left of Rectangle 3
        x[2] + w[2] + r <= x[3] # Rectangle 2 is at the left of Rectangle 3
        x[3] + w[3] + r <= x[5] # Rectangle 3 is at the left of Rectangle 5
        x[4] + w[4] + r <= x[5] # Rectangle 4 is at the left of Rectangle 5
        y[4] + h[4] + r <= H    # No rectangles on top of Rectangle 4
        y[5] + h[5] + r <= H    # No rectangles on top of Rectangle 5
        y[2] + h[2] + r <= y[1] # Rectangle 2 is below Rectangle 1
        y[1] + h[1] + r <= y[4] # Rectangle 1 is below Rectangle 4
        y[3] + h[3] + r <= y[4] # Rectangle 3 is below Rectangle 4
        w .<= 5*h               # Aspect ratio constraint
        h .<= 5*w               # Aspect ratio constraint
        A .<= h .*  w           # Area constraint
    end)

    @objective(floor_planning, Min, W + H)

    optimize!(floor_planning)

    @show objective_value(floor_planning);

    D = DataFrame(x = value.(x), y = value.(y), x2 = value.(x) .+ value.(w), y2 = value.(y) .+ value.(h))
    plt = plot(D, xmin = :x, ymin = :y, xmax = :x2, ymax = :y2, Geom.rect)
    push!(figs, plt)
end

#+

draw(SVG(6inch, 6inch), vstack(hstack(figs[1], figs[2]), hstack(figs[3], figs[4])))

#' ### References
#'
#' 1. Boyd, S., & Vandenberghe, L. (2004). Convex Optimization. Cambridge: Cambridge University Press. doi:10.1017/CBO9780511804441


text:
using Test
using CartesianMesh
using PolynomialBasis
# using Revise
using CutCellDG
include("../useful_routines.jl")

basis = LagrangeTensorProductBasis(2,1)
mesh = UniformMesh([0.,0.],[2.,1.],[2,2])
cellmaps = CutCellDG.construct_cell_maps(mesh)

points = interpolation_points(basis)
nodalcoordinates = CutCellDG.dg_nodal_coordinates(cellmaps,points)
testnodalcoordinates = [0.  0.  1.  1.  0.  0.  1.  1.  1.  1.  2.  2.  1.  1.  2.  2.
                        0.  .5  0.  .5  .5  1.  .5  1.  0.  .5  0.  .5  .5  1.  .5  1.]
@test allapprox(nodalcoordinates,testnodalcoordinates)

nodalconnectivity = CutCellDG.dg_nodal_connectivity(4,4)
testnodalconnectivity = [1  5  9   13
                         2  6  10  14
                         3  7  11  15
                         4  8  12  16]
@test allequal(nodalconnectivity,testnodalconnectivity)

basis = LagrangeTensorProductBasis(2,2)
points = interpolation_points(basis)
dgmesh = CutCellDG.DGMesh([0.,0.],[2.,1.],[2,2],points)
nodalcoordinates = CutCellDG.nodal_coordinates(dgmesh)
tn1 = [0.  0.   0.  0.5  0.5  0.5  1.  1.   1.  0.  0.   0.  0.5  0.5  0.5  1.  1.  1.
       0.  0.25 0.5 0.   0.25 0.5  0.  0.25 0.5 0.5 0.75 1.  0.5  0.75 1.   0.5 0.75 1.]
tn2 = copy(tn1)
tn2[1,:] .+= 1.0
testnodalcoordinates = hcat(tn1,tn2)
@test allapprox(nodalcoordinates,testnodalcoordinates)

nodalconnectivity = hcat([CutCellDG.nodal_connectivity(dgmesh,i) for i = 1:4]...)
testnodalconnectivity = reshape(1:36,9,:)
@test allequal(nodalconnectivity,testnodalconnectivity)


text:
#-------------------------------------------------------------------------------
# This file contains utility functions for converting undecorated source
# strings into Julia values.  For example, string->number, string unescaping, etc.

"""
Convert a Julia source code string into a number.
"""
function julia_string_to_number(str::AbstractString, kind)
    str = replace(replace(str, '_'=>""), '−'=>'-')
    if kind == K"Integer"
        x = Base.tryparse(Int, str)
        if Int === Int32 && isnothing(x)
            x = Base.tryparse(Int64, str)
        end
        if isnothing(x)
            x = Base.tryparse(Int128, str)
            if isnothing(x)
                x = Base.parse(BigInt, str)
            end
        end
        return x
    elseif kind == K"Float"
        if !startswith(str,"0x") && 'f' in str
            # This is kind of awful. Should we have a separate Float32 literal
            # type produced by the lexer?  The `f` suffix is nonstandard after all.
            return Base.parse(Float32, replace(str, 'f'=>'e'))
        else
            return Base.parse(Float64, str)
        end
    elseif kind == K"HexInt"
        ndigits = length(str)-2
        return ndigits <= 2  ? Base.parse(UInt8, str)   :
               ndigits <= 4  ? Base.parse(UInt16, str)  :
               ndigits <= 8  ? Base.parse(UInt32, str)  :
               ndigits <= 16 ? Base.parse(UInt64, str)  :
               ndigits <= 32 ? Base.parse(UInt128, str) :
               Base.parse(BigInt, str)
    elseif kind == K"BinInt"
        ndigits = length(str)-2
        return ndigits <= 8   ? Base.parse(UInt8, str)   :
               ndigits <= 16  ? Base.parse(UInt16, str)  :
               ndigits <= 32  ? Base.parse(UInt32, str)  :
               ndigits <= 64  ? Base.parse(UInt64, str)  :
               ndigits <= 128 ? Base.parse(UInt128, str) :
               Base.parse(BigInt, str)
    elseif kind == K"OctInt"
        ndigits = length(str)-2
        x = Base.tryparse(UInt64, str)
        if isnothing(x)
            x = Base.tryparse(UInt128, str)
            if isnothing(x)
                x = Base.parse(BigInt, str)
            elseif ndigits > 43
                x = BigInt(x)
            end
        else
            x = ndigits <= 3  && x <= typemax(UInt8)  ? UInt8(x)   :
                ndigits <= 6  && x <= typemax(UInt16) ? UInt16(x)  :
                ndigits <= 11 && x <= typemax(UInt32) ? UInt32(x)  :
                ndigits <= 22                         ? x          :
                ndigits <= 43                         ? UInt128(x) :
                BigInt(x)
        end
        return x
    end
end

#-------------------------------------------------------------------------------
is_indentation(c) = c == ' ' || c == '\t'

"""
Process Julia source code escape sequences for raw strings
"""
function unescape_raw_string(io::IO, str::AbstractString, is_cmd::Bool, dedent::Integer, skip_initial_newline::Bool)
    delim = is_cmd ? '`' : '"'
    i = firstindex(str)
    lastidx = lastindex(str)
    while i <= lastidx
        c = str[i]
        if c != '\\'
            if c == '\r'
                # convert literal \r and \r\n in strings to \n (issue #11988)
                if i+1 <= lastidx && str[i+1] == '\n'
                    i += 1
                end
                c = '\n'
            end
            if c == '\n'
                if i > 1 || !skip_initial_newline
                    write(io, c)
                end
                if i+1 <= lastidx && str[i+1] != '\n' && str[i+1] != '\r'
                    i += dedent
                end
            else
                write(io, c)
            end
            i = nextind(str, i)
            continue
        end
        # Process \ escape sequences
        j = i
        while j <= lastidx && str[j] == '\\'
            j += 1
        end
        nbackslash = j - i
        if (j <= lastidx && str[j] == delim) || j > lastidx
            # Backslashes before a delimiter must also be escaped
            nbackslash = div(nbackslash,2)
        end
        for k = 1:nbackslash
            write(io, '\\')
        end
        i = j
        if i <= lastidx
            write(io, str[i])
            i = nextind(str, i)
        end
    end
end

"""
Process Julia source code escape sequences for non-raw strings.
`str` should be passed without delimiting quotes.
"""
function unescape_julia_string(io::IO, str::AbstractString, dedent::Integer, skip_initial_newline::Bool)
    i = firstindex(str)
    lastidx = lastindex(str)
    while i <= lastidx
        c = str[i]
        if c != '\\'
            if c == '\r'
                # convert literal \r and \r\n in strings to \n (issue #11988)
                if i+1 <= lastidx && str[i+1] == '\n'
                    i += 1
                end
                c = '\n'
            end
            if c == '\n'
                if i > 1 || !skip_initial_newline
                    write(io, c)
                end
                if i+1 <= lastidx && str[i+1] != '\n' && str[i+1] != '\r'
                    i += dedent
                end
            else
                write(io, c)
            end
            i = nextind(str, i)
            continue
        end
        # Process \ escape sequences.  See also Base.unescape_string which some
        # of this code derives from (but which disallows \` \' \$)
        i += 1
        if i > lastidx
            break
        end
        c = str[i]
        if c == 'x' || c == 'u' || c == 'U'
            n = k = 0
            m = c == 'x' ? 2 :
                c == 'u' ? 4 : 8
            while (k += 1) <= m && i+1 <= lastidx
                nc = str[i+1]
                n = '0' <= nc <= '9' ? n<<4 + (nc-'0') :
                    'a' <= nc <= 'f' ? n<<4 + (nc-'a'+10) :
                    'A' <= nc <= 'F' ? n<<4 + (nc-'A'+10) : break
                i += 1
            end
            if k == 1 || n > 0x10ffff
                u = m == 4 ? 'u' : 'U'
                throw(ArgumentError("invalid $(m == 2 ? "hex (\\x)" :
                                    "unicode (\\$u)") escape sequence"))
            end
            if m == 2 # \x escape sequence
                write(io, UInt8(n))
            else
                print(io, Char(n))
            end
        elseif '0' <= c <= '7'
            k = 1
            n = c-'0'
            while (k += 1) <= 3 && i+1 <= lastidx
                c = str[i+1]
                n = ('0' <= c <= '7') ? n<<3 + c-'0' : break
                i += 1
            end
            if n > 255
                throw(ArgumentError("octal escape sequence out of range"))
            end
            write(io, UInt8(n))
        elseif c == '\n' || c == '\r'
            # Remove \n \r and \r\n newlines + indentation following \
            if c == '\r' && i < lastidx && str[i+1] == '\n'
                i += 1
            end
            while i < lastidx && is_indentation(str[i+1])
                i += 1
            end
        else
            u = # C escapes
                c == 'n' ? '\n' :
                c == 't' ? '\t' :
                c == 'r' ? '\r' :
                c == 'e' ? '\e' :
                c == 'b' ? '\b' :
                c == 'f' ? '\f' :
                c == 'v' ? '\v' :
                c == 'a' ? '\a' :
                # Literal escapes allowed in Julia source
                c == '\\' ? '\\' :
                c == '\'' ? '\'' :
                c == '"' ? '"' :
                c == '$' ? '$' :
                c == '`' ? '`' :
                throw(ArgumentError("Invalid escape sequence \\$c"))
            write(io, u)
        end
        i = nextind(str, i)
    end
end

function unescape_julia_string(str::AbstractString, is_cmd::Bool,
                               is_raw::Bool, dedent::Integer=0,
                               skip_initial_newline=false)
    io = IOBuffer()
    if is_raw
        unescape_raw_string(io, str, is_cmd, dedent, skip_initial_newline)
    else
        unescape_julia_string(io, str, dedent, skip_initial_newline)
    end
    String(take!(io))
end

# Compute length of longest common prefix of mixed spaces and tabs, in
# characters (/bytes).
#
# Initial whitespace is never regarded as indentation in any triple quoted
# string chunk, as it's always preceded in the source code by a visible token
# of some kind; either a """ delimiter or $() interpolation.
#
# This pass runs *before* normalization of newlines so that
# unescaping/normalization can happen in a single pass.
#
# TODO: Should we do triplequoted string splitting as part of the main parser?
# It would be conceptually clean if the trivial whitespace was emitted as
# syntax trivia.
#
# flisp: triplequoted-string-indentation-
function triplequoted_string_indentation(strs, is_raw)
    if isempty(strs)
        return 0
    end
    refstr = SubString(strs[1], 1, 0)
    reflen = -1
    for str in strs
        i = 1
        lastidx = lastindex(str)
        while i <= lastidx
            c = str[i]
            if c == '\\' && !is_raw
                # Escaped newlines stop indentation detection for the current
                # line but do not start detection of indentation on the next
                # line
                if i+1 <= lastidx
                    if str[i+1] == '\n'
                        i += 1
                    elseif str[i+1] == '\r'
                        i += 1
                        if i+1 <= lastidx && str[i+1] == '\n'
                            i += 1
                        end
                    end
                end
            elseif c == '\n' || c == '\r'
                while i <= lastidx
                    c = str[i]
                    (c == '\n' || c == '\r') || break
                    i += 1
                end
                if i <= lastidx
                    # At this point we've found the start of a nonempty line.
                    if reflen < 0
                        # Find indentation we'll use as a reference
                        j = i-1
                        while j < lastidx && is_indentation(str[j+1])
                            j += 1
                        end
                        refstr = SubString(str, i, j)
                        reflen = j - i + 1
                        if j > i
                            i = j
                        end
                    else
                        # Matching indentation with reference, shortening
                        # length if necessary.
                        j = i-1
                        while j+1 <= lastidx && j-i+2 <= reflen
                            if str[j+1] != refstr[j-i+2]
                                break
                            end
                            j += 1
                        end
                        if j-i+1 < reflen
                            reflen = j-i+1
                        end
                        if j > i
                            i = j
                        end
                    end
                else
                    # A newline directly before the end of the string means a
                    # delimiter was in column zero, implying zero indentation.
                    reflen = 0
                end
            end
            i <= lastidx || break
            i = nextind(str, i)
        end
    end
    max(reflen, 0)
end

function process_triple_strings!(strs, is_raw)
    if isempty(strs)
        return strs
    end
    dedent = triplequoted_string_indentation(strs, is_raw)
    for i = 1:length(strs)
        strs[i] = unescape_julia_string(strs[i], false, is_raw, dedent, i==1)
    end
    strs
end

#-------------------------------------------------------------------------------
# Unicode normalization. As of Julia 1.8, this is part of Base and the Unicode
# stdlib under the name `Unicode.julia_chartransform`. See
# https://github.com/JuliaLang/julia/pull/42561
#
# To allow use on older Julia versions, we reproduce that logic here.

# static wrapper around user callback function
utf8proc_custom_func(codepoint::UInt32, callback::Any) =
    UInt32(callback(codepoint))::UInt32

function utf8proc_decompose(str, options, buffer, nwords, chartransform::T) where T
    ret = ccall(:utf8proc_decompose_custom, Int, (Ptr{UInt8}, Int, Ptr{UInt8}, Int, Cint, Ptr{Cvoid}, Ref{T}),
                str, sizeof(str), buffer, nwords, options,
                @cfunction(utf8proc_custom_func, UInt32, (UInt32, Ref{T})), chartransform)
    ret < 0 && utf8proc_error(ret)
    return ret
end

function utf8proc_map(str::Union{String,SubString{String}}, options::Integer, chartransform=identity)
    nwords = utf8proc_decompose(str, options, C_NULL, 0, chartransform)
    buffer = Base.StringVector(nwords*4)
    nwords = utf8proc_decompose(str, options, buffer, nwords, chartransform)
    nbytes = ccall(:utf8proc_reencode, Int, (Ptr{UInt8}, Int, Cint), buffer, nwords, options)
    nbytes < 0 && utf8proc_error(nbytes)
    return String(resize!(buffer, nbytes))
end

const _julia_charmap = Dict{UInt32,UInt32}(
    0x025B => 0x03B5,
    0x00B5 => 0x03BC,
    0x00B7 => 0x22C5,
    0x0387 => 0x22C5,
    0x2212 => 0x002D,
)

julia_chartransform(codepoint::UInt32) = get(_julia_charmap, codepoint, codepoint)

function normalize_identifier(str)
    flags = Base.Unicode.UTF8PROC_STABLE | Base.Unicode.UTF8PROC_COMPOSE
    utf8proc_map(str, flags, julia_chartransform)
end


text:
using ResultTypes
using Test

struct BarError <: Exception end
struct FooError <: Exception end

@testset "ResultTypes" begin

@testset "Result" begin
    @testset "Basic result" begin
        x = Result(2)
        @test unwrap(x) === 2
        @test_throws ErrorException unwrap_error(x)
        @test ResultTypes.iserror(x) === false
        @test x isa Result{Int, Exception}

        # on unwrapped already
        y = unwrap(x)
        @test unwrap(y) === y
    end

    @testset "Result with error type" begin
        x = Result(2, DivideError)
        @test unwrap(x) === 2
        @test_throws ErrorException unwrap_error(x)
        @test ResultTypes.iserror(x) === false
        @test x isa Result{Int, DivideError}
    end

    @testset "Malformed result" begin
        x = Result{Int, Exception}(nothing, nothing)
        @test_throws ErrorException unwrap(x)
        @test_throws ErrorException unwrap_error(x)
    end
end

@testset "ErrorResult" begin
    @testset "Basic error" begin
        x = ErrorResult(Int, "Basic Error")
        @test_throws ErrorException unwrap(x)
        e = unwrap_error(x)
        @test isa(e, ErrorException)
        @test e.msg == "Basic Error"
        @test ResultTypes.iserror(x) === true

        # directly on exceptions
        @test ResultTypes.iserror(e)
        @test unwrap_error(e) === e
        @test !ResultTypes.iserror(2)
    end

    @testset "Empty error" begin
        x = ErrorResult(Int)
        @test_throws ErrorException unwrap(x)
        e = unwrap_error(x)
        @test isa(e, ErrorException)
        @test ResultTypes.iserror(x) === true
        @test x isa Result{Int, ErrorException}

        x = ErrorResult()
        @test_throws ErrorException unwrap(x)
        e = unwrap_error(x)
        @test isa(e, ErrorException)
        @test ResultTypes.iserror(x) === true
        @test x isa Result{Any, ErrorException}
    end

    @testset "Typed Error" begin
        x = ErrorResult(Int, DivideError())
        @test_throws DivideError unwrap(x)
        @test isa(unwrap_error(x), DivideError)
        @test ResultTypes.iserror(x) === true
        @test x isa Result{Int, DivideError}
    end
end

@testset "Convert" begin
    @testset "From Result" begin
        x = Result(2)
        @test unwrap(Int, x) === 2
        @test unwrap(Float64, x) === 2.0
        @test_throws MethodError unwrap(String, x)
    end

    @testset "From ErrorResult" begin
        x = ErrorResult(Int, "Foo")
        @test_throws ErrorException unwrap(Int, x)
    end

    @testset "To Result" begin
        x = convert(Result{Int, ErrorException}, 2.0)
        @test unwrap(x) === 2
        @test isa(x, Result{Int, ErrorException})
    end

    @testset "To ErrorResult" begin
        x = convert(Result{Int, DivideError}, DivideError())
        @test unwrap_error(x) == DivideError()
        @test isa(x, Result{Int, DivideError})
    end

    @testset "Result to Result" begin
        r = Result{Int, DivideError}(DivideError())
        x = convert(Result{Int, DivideError}, r)
        @test r == x
        @test r === x

        r = Result{Int, ErrorException}(2)
        x = convert(Result{Int, ErrorException}, r)
        @test r == x

        r = Result{Int, KeyError}(KeyError("foo"))
        x = convert(Result{Real, Exception}, r)
        @test x isa Result{Real, Exception}
        @test_throws KeyError unwrap(x)

        r = Result{Int, KeyError}(2)
        x = convert(Result{Real, Exception}, r)
        @test x isa Result{Real, Exception}
        @test unwrap(x) == 2
    end

    @testset "Promote" begin
        @test eltype([Result(2), ErrorResult(Int, DivideError())]) == Result{Int, Exception}
        @test eltype([Result(2), ErrorResult(DivideError())]) == Result{Any, Exception}
        @test eltype([Result(2.0), ErrorResult(Int, KeyError(:foo))]) == Result{Float64, Exception}
        @test eltype([Result(2.0, KeyError), ErrorResult(Int, KeyError(:foo))]) == Result{Float64, KeyError}

        # only for coverage, this is called above
        @test promote_rule(Result{Float64, Exception}, Result{Int, KeyError}) == Result{Float64, Exception}
    end

    @testset "Ambiguities" begin
        ambs = filter(Test.detect_ambiguities(Base, ResultTypes)) do amb
            undeprecated = all(amb) do m
                # method is deprecated, here or in Base
                m.file != Symbol("deprecated.jl")
            end

            ours = any(amb) do m
                m.module == ResultTypes
            end

            undeprecated && ours
        end

        @test isempty(ambs)
    end
end

@testset "Return" begin
    let
        function foo(x::Int, y::Int)::Result{Int, DivideError}
            try
                return div(x, y)
            catch e
                return e
            end
        end

        x = foo(3, 2)
        @test isa(x, Result{Int, DivideError})
        @test unwrap(x) === 1

        y = foo(1, 0)
        @test isa(y, Result{Int, DivideError})
        @test unwrap_error(y) == DivideError()
    end
end

@testset "@try" begin
    Base.convert(::Type{FooError}, ::BarError) = FooError()

    function isbar(x)::Result{Bool, BarError}
        x == "foo" && return BarError()
        return x == "bar"
    end

    function foo(x)::Result{Int, FooError}
        bar = @try isbar(x)
        return bar ? 42 : 13
    end

    foobar = foo("bar")
    @test isa(foobar, Result{Int, FooError})
    @test unwrap(foobar) === 42

    foofoo = foo("foo")
    @test isa(foofoo, Result{Int, FooError})
    @test unwrap_error(foofoo) === FooError()

    football = foo("tball")
    @test isa(football, Result{Int, FooError})
    @test unwrap(football) === 13

    function foo2(x)::Result{Int, ErrorException}
        bar = @try(isbar(x), ErrorException("Got a BarError"))
        return bar ? 43 : 14
    end

    foobar = foo2("bar")
    @test isa(foobar, Result{Int, ErrorException})
    @test unwrap(foobar) === 43

    foofoo = foo2("foo")
    @test isa(foofoo, Result{Int, ErrorException})
    @test unwrap_error(foofoo) === ErrorException("Got a BarError")

    football = foo2("tball")
    @test isa(football, Result{Int, ErrorException})
    @test unwrap(football) === 14
end

@testset "Show" begin
    @testset "Result" begin
        x = Result(2)

        @test isa(string(x), String)
    end

    @testset "ErrorResult" begin
        x = ErrorResult(Int, "Basic Error")
        @test isa(string(x), String)
    end
end

end


text:
import XLSX

#
# examples from docstrings
#

data_directory = joinpath(dirname(@__FILE__), "..", "data")

v = XLSX.readdata(joinpath(data_directory, "myfile.xlsx"), "mysheet", "A1:B4")

f = XLSX.openxlsx(joinpath(data_directory, "myfile.xlsx"))
sheet = f["mysheet"]
v = sheet["A1:B4"]
close(f)

xf = XLSX.openxlsx(joinpath(data_directory, "myfile.xlsx"))
sheet = xf["mysheet"]
cell = XLSX.getcell(sheet, "A1")
close(xf)

xf = XLSX.openxlsx(joinpath(data_directory, "myfile.xlsx"))
sheet = xf["mysheet"]
left = 1  # 1st column
right = 4 # 4th column
for sheetrow in XLSX.eachrow(sheet)
    for column in left:right
        cell = XLSX.getcell(sheetrow, column)

        # do something with cell
    end
end

using DataFrames, XLSX
xf = XLSX.openxlsx(joinpath(data_directory, "myfile.xlsx"))
df = DataFrame(XLSX.gettable(xf["mysheet"])...)
close(xf)

XLSX.decode_column_number("D")
XLSX.encode_column_number(4)

cn = XLSX.CellRef("AB1")
XLSX.row_number(cn); # will print 1
XLSX.column_number(cn); # will print 28
string(cn); # will print out AB1

cn = XLSX.ref"AB1"
XLSX.row_number(cn); # will print 1
XLSX.column_number(cn); # will print 28
string(cn); # will print out AB1
cr = XLSX.range"A1:C4"

#
# examples from README.md
#

xf = XLSX.openxlsx(joinpath(data_directory, "myfile.xlsx"))
XLSX.sheetnames(xf)
sh = xf["mysheet"]
sh["B2"] # access a cell value
sh["A2:B4"] # access a range
XLSX.readdata(joinpath(data_directory, "myfile.xlsx"), "mysheet", "A2:B4") # shorthand for all above
xf["mysheet!A2:B4"] # you can also query values from a file reference
xf["NAMED_CELL"] # you can even read named ranges
xf["mysheet!A:B"] # Column ranges are also supported
sh[:] # all data inside worksheet's dimension
XLSX.getdata(sh) # same as sh[:]
close(xf)

using DataFrames, XLSX

df = DataFrame(XLSX.readtable(joinpath(data_directory, "myfile.xlsx"), "mysheet")...)

using DataFrames, XLSX

df = DataFrame(XLSX.readtable(joinpath(data_directory, "myfile.xlsx"), "mysheet")...)

f = XLSX.openxlsx(joinpath(data_directory, "myfile.xlsx"), enable_cache=false)
sheet = f["mysheet"]
for r in XLSX.eachrow(sheet)
    # r is a `SheetRow`, values are read using column references
    rn = XLSX.row_number(r) # `SheetRow` row number
    v1 = r[1]    # will read value at column 1
    v2 = r["B"]  # will read value at column 2
end

for r in XLSX.eachtablerow(sheet)
	# r is a `TableRow`, values are read using column labels or numbers
	rn = XLSX.row_number(r) # `TableRow` row number
	v1 = r[1] # will read value at table column 1
	v2 = r[:HeaderB] # will read value at column labeled `:HeaderB`
end

import DataFrames, XLSX
df = DataFrames.DataFrame(integers=[1, 2, 3, 4], strings=["Hey", "You", "Out", "There"], floats=[10.2, 20.3, 30.4, 40.5])
XLSX.writetable("df.xlsx", DataFrames.columns(df), DataFrames.names(df))

df1 = DataFrames.DataFrame(COL1=[10,20,30], COL2=["Fist", "Sec", "Third"])
df2 = DataFrames.DataFrame(AA=["aa", "bb"], AB=[10.1, 10.2])

XLSX.writetable("report.xlsx", REPORT_A=( DataFrames.columns(df1), DataFrames.names(df1) ), REPORT_B=( DataFrames.columns(df2), DataFrames.names(df2) ))


text:
"""
    WeightedParticleRepresentation{S} <: ParticleRepresentation{S} <: AbstractFilterRepresentation{S}

Abstract type for the representation of the conditional distribution over the hidden state by weighted particles (samples).
"""
abstract type WeightedParticleRepresentation{S} <: ParticleRepresentation{S} end

"""
    get_weight(ensemble::WeightedParticleRepresentation, i)

Return the weight of the particle with index `i` from `ensemble`.
"""
function get_weight(ensemble::WeightedParticleRepresentation, i) end

"""
    list_of_weights(ensemble::WeightedParticleRepresentation)

Return a list of particle weights from `ensemble`.
"""
list_of_weights(ens::WeightedParticleRepresentation) = [get_weight(ens, i) for i in 1:no_of_particles(ens)]

"""
    sum_of_weights(ensemble::WeightedParticleRepresentation)

Return the sum of the importance weights in `ensemble`.
"""
function sum_of_weights(ens::WeightedParticleRepresentation)
    sum = 0.
    for i in 1:no_of_particles(ens)
        sum += get_weight(ens, i)
    end
    return sum
end

"""
    resample!(ensemble::WeightedParticleRepresentation)

Resample the positions by drawing a weighted sample with replacement.
"""
function resample!(ens::WeightedParticleRepresentation) end

function eff_no_of_particles(ensemble::WeightedParticleRepresentation)
    sum = 0.
    @inbounds for i in 1:no_of_particles(ensemble)
        sum += get_weight(ensemble, i).^2
    end
    return sum_of_weights(ensemble)^2/sum
end

dim(ensemble::WeightedParticleRepresentation) = (particle_dim(ensemble) + 1) * no_of_particles(ensemble)

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
# # Raw data plots
#
# ### Summary of instrumentation
#
# For the paper we only used the measurements of two DCX-22-CTD sensors, 205309-100mH2O and 207265-300mH2O (referred to as CTD-309 and CTD-265, respectively).
# We installed a third one, 205145-10mH2O (CTD-145), on 20 and 21 of August but could not use it for further calculations (see `derived_quantities.pdf`).
# Additionally, we monitored the pressure in the artificial moulin over night and the pressure in the supplying water stream and the atmosphere.
# From 8 August to 19 August the instruments were installed in AM15 and on 20 and 21 August in AM13.
#
# - 3x DCX-22-CTD: CTD-265, CTD-309 and CTD-145 (the latter was only operating on August 20 and 21), installed at the following depths:
#
#
#   | Day         | 08.08. | 09.-12.08. | 13.-14.08., 17.-18.08. | 19.08. | 20.-21.08. |
#   |------------ | ------ |------------| ---------------------- | ------ | ---------- |
#   | CTD-265 [m] |    160 |        167 |                    167 |    167 |        185 |
#   | CTD-309 [m] |     60 |        117 |                    127 |    137 |         85 |
#   | CTD-145 [m] |      - |          - |                      - |      - |         10 |
#
#
# - 1x DCX-22 pressure sensor (over night), installed at:
#
#   | Day       | 07.08. | 08.08. | 09.08. | 11.-13.08., 17.-18.08. | 19.08. | 20.08. |
#   |---------- | ------ |--------| ------ | ---------------------- | ------ | ------ |
#   | Depth [m] |     90 |     75 |    140 |                    150 |    167 |      ? |
#
#
# - 1x DCX-22AA stage sensor (installed a few metres upstream of artificial moulins), measuring pressure and temperature in the stream and the atmosphere

using Rhone2020  #hide
const R = Rhone2020  #hide

using PyPlot #hide
include(joinpath(@__DIR__, "../read_measurements.jl")); #hide

# ### Example traces of conductivity
#
# These traces were manually cut out from the conductivity time series and used to calculate the discharges and the mean velocity.
#
# This example shows the traces from the fifth injection on August 21.

R.plot_example_traces((ctd145, ctd309, ctd265), indices, "2108", 5) #hide

# ### Time series of individual days
#
# Conductivity, pressure and temperature data from individual days. Note that where the conductivity is greater than the plot limit the sensor was supersaturated.
# Calibrations were carried out August 5, 10, 13 and 20.
# The temperature data are already corrected according to the zero-degree calibration.
# Shaded areas in the temperature plots show the uncertainty (+/- 0.05°C).

R.plot_CTD("0808", ctd309, ctd265, ctd145, e_T, 15500:28640) #hide
#+
R.plot_CTD("0908", ctd309, ctd265, ctd145, e_T, 8000:28640) #hide
#+
R.plot_CTD("1008", ctd309, ctd265, ctd145, e_T, 11500:20500) #hide
#+
R.plot_CTD("1108", ctd309, ctd265, ctd145, e_T, 4400:25700) #hide
#+
R.plot_CTD("1208", ctd309, ctd265, ctd145, e_T, 3400:26000) #hide
#+
R.plot_CTD("1308", ctd309, ctd265, ctd145, e_T, 4900:15900) #hide
#+
R.plot_CTD("1408", ctd309, ctd265, ctd145, e_T, 4100:20000) #hide
#+
R.plot_CTD("1708", ctd309, ctd265, ctd145, e_T, 9400:24000) #hide
#+
R.plot_CTD("1808", ctd309, ctd265, ctd145, e_T, 1700:28640) #hide
#+
R.plot_CTD("1908", ctd309, ctd265, ctd145, e_T, 1100:28640) #hide
#+
R.plot_CTD("2008", ctd309, ctd265, ctd145, e_T, 3600:25800) #hide
#+
R.plot_CTD("2108", ctd309, ctd265, ctd145, e_T, 1:28640) #hide


text:
# Copyright 2020-2021 John T. Foster
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
using Test
using TaD
import TaD: evaluate

@testset "evaluate bspline basis functions" begin
    kv = [0, 0, 0, 1, 2, 3, 4, 4, 5, 5, 5.0]
    @test evaluate(BSplineBasis(kv, 2), 5 / 2)[1] ≈ 1 / 8
    @test evaluate(BSplineBasis(kv, 2), 5 / 2)[2] ≈ 3 / 4
    @test evaluate(BSplineBasis(kv, 2), 5 / 2)[3] ≈ 1 / 8
end
@testset "evaluate bspline basis function derivatives" begin
    kv = [0, 0, 0, 0, 0.5, 1, 1, 1, 1]
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[1, 1] ≈ 0.25
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[1, 2] ≈ 0.5
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[1, 3] ≈ 0.25
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[1, 4] ≈ 0.0
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[2, 1] ≈ -1.5
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[2, 2] ≈ 0.0
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[2, 3] ≈ 1.5
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[2, 4] ≈ 0.0
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[3, 1] ≈ 6.0
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[3, 2] ≈ -12.0
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[3, 3] ≈ 6.0
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[3, 4] ≈ 0.0
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[4, 1] ≈ -12.0
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[4, 2] ≈ 48.0
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[4, 3] ≈ -84.0
    @test evaluate(BSplineBasis(kv, 3, k=3), 0.5)[4, 4] ≈ 48.0
end
@testset "evaluate bspline curve and derivatives" begin
    kv = [0, 0, 0, 1, 1, 1]
    b = BSplineBasis(kv, 2, k=1)
    cp = [0 0; 0.5 1; 1 0]
    c = BSplineCurve(b, cp)
    @test evaluate(c, 0.05)[1, 1] ≈ 0.05
    @test evaluate(c, 0.05)[1, 2] ≈ 0.095
    @test evaluate(c, 0.05)[2, 1] ≈ 1.0
    @test evaluate(c, 0.05)[2, 2] ≈ 1.8
    @test evaluate(c, 0.5)[1, 1] ≈ 0.5
    @test evaluate(c, 0.5)[1, 2] ≈ 0.5
    @test evaluate(c, 0.5)[2, 1] ≈ 1.0
    @test evaluate(c, 0.5)[2, 2] ≈ 0.0
end


text:
function default_encoding()
  grn_encoding = ccall((:grn_get_default_encoding, libgroonga),
                       Ptr{Void}, ())
  encoding = ccall((:grn_encoding_to_string, libgroonga),
                   Ptr{UInt8}, (Ptr{Void},), grn_encoding)
  bytestring(encoding)
end


text:
"""
`PeakResults`

This struct holds peak properties.

**Fields**
- `indices`: Vector of peak indices.
- `peaks`: Vector of peak values.
- `locations`: Vector of peak locations.
- `prominences`: Vector of peak prominences.
- `width_bounds`: Matrix with two columns indicating the beginnings and the ends
    of the peak widths.
"""
struct PeakResults
    indices::Vector  # Peak indices
    peaks::Vector  # Peak values
    locations::Vector  # Peak locations
    prominences::Vector  # Prominences
    width_bounds::Matrix  # Width bounds
end

"""
`peakheights(pr::PeakResults)`

Return peak values (local maxima).
"""
function peakheights(pr::PeakResults)
    return pr.peaks
end

"""
`peaklocations(pr::PeakResults)`

Return peak locations.
"""
function peaklocations(pr::PeakResults)
    return pr.locations
end

"""
`peakprominences(pr::PeakResults)`

Return peak prominences.
"""
function peakprominences(pr::PeakResults)
    return pr.prominences
end

"""
`peakwidths(pr::PeakResults)`

Return peak widths.
"""
function peakwidths(pr::PeakResults)
    return [diff(pr.width_bounds, dims = 2)...]
end

"""
`isempty(pr::PeakResults)`

Return `true` if `pr` does not contain any peaks.
"""
function isempty(pr::PeakResults)
    return isempty(pr.indices)
end

"""
`length(pr::PeakResults)`

Return the number of peaks in `pr`.
"""
function length(pr::PeakResults)
    return length(pr.indices)
end

"""
`getindex(pr::PeakResults, inds)`

Return a subset of `pr` as specified by `inds`.
"""
function getindex(pr::PeakResults, inds)
    indices = pr.indices[inds]
    peaks = pr.peaks[inds]
    locations = pr.locations[inds]
    prominences = pr.prominences[inds]
    width_bounds = pr.width_bounds[inds, :]
    return PeakResults(indices, peaks, locations, prominences, width_bounds)
end

"""
`sort(pr::PeakResults; rev::Bool=false, ref::String="height")`

Return sorted `pr` with peaks sorted from lowest to highest. Use `rev = true` to
    reverse the sorting order. Use `ref` to specify the property that is used to
    sort the peaks. The possible options are "height" (default), "prominence",
    and "location".
"""
function sort(pr::PeakResults; rev::Bool=false, ref::String="height")

    if ref == "height"
        peaks = pr.peaks
        inds = sortperm(peaks; rev=rev)
    elseif ref == "prominence"
        prominences = pr.prominences
        inds = sortperm(prominences; rev=rev)
    elseif ref == "location"
        locations = pr.locations
        inds = sortperm(locations; rev=rev)
    else
        throw(ArgumentError("Unknown value of ref: $(ref)"))
    end
    pr_sorted = pr[inds]
    return pr_sorted

end

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
using Printf
using Oceananigans
using Oceananigans.BoundaryConditions
using Oceananigans.Fields
using Oceananigans.OutputWriters
using Oceananigans.Diagnostics
using Oceananigans.Utils
using Oceananigans.AbstractOperations
using Vizinanigans

## Create Oceananigans Model
arch = CPU()
FT   = Float64

end_time = 100*365day
const Lx = 1000kilometer #
const Ly = 1000kilometer #
const Lz = 3kilometer    #

# Rough resolution
Δx = Δy = scale * 250meter # 5km
Δz = 100meter

const Nx = round(Int, Lx/ Δx / 16) * 16
const Ny = round(Int, Ly/ Δy / 16) * 16
const Nz = round(Int, Lz/ Δz / 16) * 16

topology = (Periodic, Bounded, Bounded)
grid = RegularCartesianGrid(topology=topology, size=(Nx, Ny, Nz), x=(0, Lx), y=(0, Ly), z=(-Lz, 0))

const f = -1e-4
const β = 1 * 10^(-11)
coriolis = FPlane(FT, f=f)
coriolis = BetaPlane(FT, f₀ = f, β = β)

α = 2e-4  # Thermal expansion coefficient [K⁻¹]
eos = LinearEquationOfState(FT, α=α, β=0)
buoyancy = BuoyancyTracer()

κh = 0.5e-5
νh = 12.0
κv = 0.5e-5
νv = 3e-4
closure = AnisotropicDiffusivity(νx = νh, νy = νh, νz =νv, κx = κh, κy = κh, κz=κv)
ΔB = 10 * 2e-3
h = 1000.0 # [m] relexaction profile scale
# Initial Conditions
U₀(x,y,z) = (z + Lz)/(f + β * y) * (ΔB / Ly)
B₀(x, y, z) = ΔB * (y / Ly +  ( exp(z/h) - exp(-Lz/h) ) / (1 - exp(-Lz/h)) - 1)

## checkpointing
model = IncompressibleModel(
            architecture = arch,
                float_type = FT,
                    grid = grid,
                coriolis = coriolis,
                buoyancy = buoyancy,
                closure = closure,
                tracers = (:b,),
)
set!(model, u = U₀, b=B₀)

## check Makie
visualize(model)

max_stars_repo_name: glwagner/Visualizanigans

max_issues_repo_name: glwagner/Visualizanigans

max_forks_repo_name: glwagner/Visualizanigans


text:
@testset "adapt_rho_gamma_parallel" begin
#test adapt_rho_gamma_parallel.jl
#test if the parallel and serial version yield the same results

#without BLAS
x=randn(100)
p=2
i=10

y=Vector{Vector{Float64}}(undef, p)
y[1]=randn(51)
y[2]=randn(100)
y_0=Vector{Vector{Float64}}(undef, p)
y_0[1]=randn(51)
y_0[2]=randn(100)
y_old=Vector{Vector{Float64}}(undef, p)
y_old[1]=randn(51)
y_old[2]=randn(100)
l_old=Vector{Vector{Float64}}(undef, p)
l_old[1]=randn(51)
l_old[2]=randn(100)
l=Vector{Vector{Float64}}(undef, p)
l[1]=randn(51)
l[2]=randn(100)
d_l=Vector{Vector{Float64}}(undef, p)
d_l[1]=randn(51)
d_l[2]=randn(100)
l_0=Vector{Vector{Float64}}(undef, p)
l_0[1]=randn(51)
l_0[2]=randn(100)
l_hat=Vector{Vector{Float64}}(undef, p)
l_hat[1]=randn(51)
l_hat[2]=randn(100)
l_hat_0=Vector{Vector{Float64}}(undef, p)
l_hat_0[1]=randn(51)
l_hat_0[2]=randn(100)
d_l_hat=Vector{Vector{Float64}}(undef, p)
d_l_hat[1]=randn(51)
d_l_hat[2]=randn(100)
d_H_hat=Vector{Vector{Float64}}(undef, p)
d_H_hat[1]=randn(51)
d_H_hat[2]=randn(100)
d_G_hat=Vector{Vector{Float64}}(undef, p)
d_G_hat[1]=randn(51)
d_G_hat[2]=randn(100)

rho=[1.234;10.23432]
gamma=Vector{Float64}(undef, p)
gamma[1]=1.0
gamma[2]=1.345

prox=Vector{Any}(undef, p)
prox[1] = x -> 1.0.*x
m=randn(100)
prox[2] = x -> prox_l2s!(x,rho[2],m)
TD_OP=Vector{SparseMatrixCSC{Float64,Int64}}(undef, p)
TD_OP[1] = sparse(I, 51, 100)*2.0
TD_OP[2] = sparse(I, 100, 100)
maxit=39
#log_PSDMM = log_type_PARSDMM(zeros(maxit,p-1),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),zeros(maxit),zeros(maxit),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),TimerOutput());
P_sub = Vector{Any}(undef, 1)
P_sub[1] = prox[1]
counter=12
x_hat=Vector{Vector{Float64}}(undef, p)
x_hat[1]=randn(51)
x_hat[2]=randn(100)
r_pri=Vector{Vector{Float64}}(undef, p)
r_pri[1]=randn(51)
r_pri[2]=randn(100)
s=Vector{Vector{Float64}}(undef, p)
s[1]=randn(51)
s[2]=randn(100)
s_0=Vector{Vector{Float64}}(undef, p)
s_0[1]=randn(51)
s_0[2]=randn(100)

#distributed quantities
s_d      = distribute(s)
s_0_d    = distribute(s_0)
r_pri_d  = distribute(r_pri)
x_hat_d  = distribute(x_hat)
prox_d   = distribute(prox)
P_sub_d  = distribute(P_sub)
TD_OP_d  = distribute(TD_OP)
rho_d    = distribute(rho)
gamma_d  = distribute(gamma)
y_d      = distribute(y)
y_0_d    = distribute(y_0)
y_old_d  = distribute(y_old)
l_d      = distribute(l)
d_l_d    = distribute(d_l)
l_old_d  = distribute(l_old)
l_0_d    = distribute(l_0)
l_hat_d  = distribute(l_hat)
d_l_hat_d= distribute(d_l_hat)
d_H_hat_d= distribute(d_H_hat)
d_G_hat_d=distribute(d_G_hat)
l_hat_0_d=distribute(l_hat_0)

#distributed computation
[ @spawnat pid  adapt_rho_gamma_parallel(gamma_d[:L],rho_d[:L],true,true,y_d[:L],y_old_d[:L],s_d[:L],s_0_d[:L],l_d[:L],l_hat_0_d[:L],
                                         l_0_d[:L],l_old_d[:L],y_0_d[:L],l_hat_d[:L],d_l_hat_d[:L],d_H_hat_d[:L],d_l_d[:L],d_G_hat_d[:L]) for pid in y_d.pids]

#serial computation
(rho,gamma,l_hat,d_l_hat,d_H_hat,d_l,d_G_hat)=adapt_rho_gamma(i,gamma,rho,true,true,y,y_old,s,s_0,l,l_hat_0,l_0,l_old,y_0,p,l_hat,d_l_hat,d_H_hat,d_l,d_G_hat);

#see if output is the same
@test rho_d==rho
@test gamma==gamma_d
@test isapprox(l_hat,l_hat_d,rtol=10*eps())
@test isapprox(d_l_hat,d_l_hat_d,rtol=10*eps())
@test isapprox(d_H_hat,d_H_hat_d,rtol=10*eps())
@test isapprox(d_l,d_l_d,rtol=10*eps())
@test isapprox(d_G_hat,d_G_hat_d,rtol=10*eps())

#test another related piece of code
for ii = 1:p
  l_hat[ii]   = l_old[ii] .+ rho[ii].* ( -s[ii] .+ y_old[ii] );
  copy!(l_hat_0[ii],l_hat[ii])
  copy!(y_0[ii],y[ii])
  copy!(s_0[ii],s[ii])
  copy!(l_0[ii],l[ii])
end

[@spawnat pid l_hat_d[:L] = l_old_d[:L] .+ rho_d[:L].* ( -s_d[:L] .+ y_old_d[:L] ) for pid in l_hat_d.pids]

[@spawnat pid copy!(l_hat_0_d[:L],l_hat_d[:L]) for pid in l_hat_d.pids]
[@spawnat pid copy!(y_0_d[:L],y_d[:L]) for pid in l_hat_d.pids]
[@spawnat pid copy!(s_0_d[:L],s_d[:L]) for pid in l_hat_d.pids]
[@spawnat pid copy!(l_0_d[:L],l_d[:L]) for pid in l_hat_d.pids]

@test isapprox(l_hat_0_d,l_hat_0,rtol=10*eps())
@test isapprox(l_hat,l_hat_d,rtol=10*eps())
@test isapprox(y_0_d,y_d,rtol=10*eps())
@test isapprox(s_0_d,s_d,rtol=10*eps())
@test isapprox(l_0_d,l_d,rtol=10*eps())

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
using LabelledArrays

export rossler_rossler_rossler_bidir_forced

@inline @inbounds function eom_rossler_rossler_rossler_bidir_forced(u, p, t)
    ω₁, ω₂, ω₃, c_xy, c_yx, c_zx, c_zy, a₁, a₂, a₃, b₁, b₂, b₃, a₃, c₁, c₂, c₃ = (p...,)
    x1, x2, x3, y1, y2, y3, z1, z2, z3 = (u...,)

    dx1 = -ω₁*(x2 + x3) + c_yx*(y1 - x1) + c_zx*(z1 - x1)
    dx2 = ω₁*x1 + a₁*x2
    dx3 = a₂ + x3*(x1 - a₃)

    dy1 = -ω₂*(y2 + y3) + c_xy*(x1 - y1) + c_zy*(z1 - y1)
    dy2 = ω₂*y1 + b₁*y2
    dy3 = b₂ + y3*(y1 - b₃)

    dz1 = -ω₂*(z2 + z3)
    dz2 = ω₂*z1 + c₁*z2
    dz3 = c₂ + z3*(z1 - c₃)

    return SVector{9}(dx1, dx2, dx3, dy1, dy2, dy3, dz1, dz2, dz3)
end

"""
    rossler_rossler_rossler_bidir_forced(; u0 = rand(9),
        ω₁ = 1.015, ω₂ = 0.985, ω₃ = 0.95,
        c_xy = 0.1, c_yx = 0.1,
        c_zx = 0.05, c_zy = 0.05,
        a₁ = 0.15, a₂ = 0.2, a₃ = 10,
        b₁ = 0.15, b₂ = 0.2, b₃ = 10,
        c₁ = 0.15, c₂ = 0.2, c₃ = 10)

Equations of motion for a system consisting of three coupled 3D Rössler systems
(``X``, ``Y``, ``Z``), giving a 9D system [1]. The external system
``Z`` influences both ``X`` and ``Y`` (controlled by `c_zx` and `c_zy`).
Simultaneously, the subsystems  ``X`` and ``Y`` bidirectionally
influences each other (controlled by `c_xy` and `c_yx`).

The ``X`` and ``Y`` subsystems are mostly synchronized for `c_xy > 0.1` or
`c_yx > 0.1`.

## Equations of motion

The dynamics is generated by the following vector field

```math
\\begin{aligned}
\\dot{x_1} &= -\\omega_1 (x_2 + x_3) + c_{yx}(y_1 - x_1) + c_{zx}(z_1 - x_1) \\\\
\\dot{x_2} &= \\omega_1 x_1 + a_1 x_2 \\\\
\\dot{x_3} &= a_2 + x_3 (x_1 - a_3) \\\\
\\dot{y_1} &= -\\omega_1 (y_2 + y_3) + c_{xy}(x_1 - y_1) + c_{zy}(z_1 - y_1) \\\\
\\dot{x_2} &= \\omega_2 y_1 + b_1 y_2 \\\\
\\dot{x_3} &= b_2 + x_3 (y_1 - b_3) \\\\
\\dot{y_1} &= -\\omega_2 (z_2  + z_3) \\\\
\\dot{x_2} &= \\omega_2 z_1 + c_1 z_2 \\\\
\\dot{x_3} &= c_2 + z_3 (z_1 - c_3)
\\end{aligned}
```

## References

1. Amigó, José M., and Yoshito Hirata. "Detecting directional couplings from
    multivariate flows by the joint distance distribution." Chaos: An
    Interdisciplinary Journal of Nonlinear Science 28.7 (2018): 075302.
"""
function rossler_rossler_rossler_bidir_forced(; u0 = rand(9),
        ω₁ = 1.015, ω₂ = 0.985, ω₃ = 0.95,
        c_xy = 0.1, c_yx = 0.1,
        c_zx = 0.05, c_zy = 0.05,
        a₁ = 0.15, a₂ = 0.2, a₃ = 10,
        b₁ = 0.15, b₂ = 0.2, b₃ = 10,
        c₁ = 0.15, c₂ = 0.2, c₃ = 10)

    p = @LArray [ω₁, ω₂, ω₃, c_xy, c_yx, c_zx, c_zy, a₁, a₂, a₃, b₁, b₂, b₃, a₃, c₁, c₂, c₃] (:ω₁, :ω₂, :ω₃, :c_xy, :c_yx, :c_zx, :c_zy, :a₁, :a₂, :a₃, :b₁, :b₂, :b₃, :a₃, :c₁, :c₂, :c₃)
    ContinuousDynamicalSystem(eom_rossler_rossler_rossler_bidir_forced, u0, p)
end

function forced_rossler_rossler_bidir_trajectory(npts;
        n_transient = 2000, dt = 0.6, sample_dt = 1,
        u0 = rand(9), ω₁ = 1.015, ω₂ = 0.985, ω₃ = 0.95,
        c_xy = 0.2, c_yx = 0.2, c_zx = 0.05, c_zy = 0.05,
        a₁ = 0.15, a₂ = 0.2, a₃ = 10,
        b₁ = 0.15, b₂ = 0.2, b₃ = 10,
        c₁ = 0.15, c₂ = 0.2, c₃ = 10)

    s = rossler_rossler_rossler_bidir_forced(u0 = u0,
        c_xy = c_xy, c_yx = c_yx,
        c_zx = c_zx, c_zy = c_zy,
        ω₁ = ω₁, ω₂ = ω₂, ω₃ = ω₃,
        a₁ = a₁, a₂ = a₂, a₃ = a₃,
        b₁ = b₁, b₂ = b₂, b₃ = b₃,
        c₁ = c₁, c₂ = c₂, c₃ = c₃)

    # the system is recorded at times t0:dt:T
    T = npts*dt*sample_dt
    o = trajectory(s, T, dt = dt, Ttr = n_transient*dt, alg = SimpleDiffEq.SimpleATsit5())[1:sample_dt:end-1, :] #alg = SimpleDiffEq.SimpleATsit5()
end

function good_rossler_rossler_rossler_bidir_forced_trajectory(npts;
        sample_dt = 1,  Ttr = 5000, dt = 0.3, # dt = 0.6 about 10 samples per period
        Da₁ = Uniform(0.12, 0.17),
        Da₂ = Uniform(0.18, 0.22),
        Da₃ = Uniform(9.0, 11.0),
        Db₁ = Uniform(0.10, 0.20),
        Db₂ = Uniform(0.18, 0.22),
        Db₃ = Uniform(9.0, 11.0),
        Dc₁ = Uniform(0.10, 0.20),
        Dc₂ = Uniform(0.18, 0.22),
        Dc₃ = Uniform(9.0, 11.0),
        Dω₁ = Uniform(0.95, 0.999), #Uniform(0.97, 1.03)
        Dω₂ = Uniform(1.001, 1.05), #Uniform(0.97, 1.03)
        Dω₃ = Uniform(0.9, 0.95), #Uniform(0.97, 1.03)

        a₁ = nothing,
        a₂ = nothing,
        a₃ = nothing,
        b₁ = nothing,
        b₂ = nothing,
        b₃ = nothing,
        c₁ = nothing,
        c₂ = nothing,
        c₃ = nothing,
        ω₁ = nothing,
        ω₂ = nothing,
        ω₃ = nothing,
        c_xy = 0.2,  c_yx = 0.2,
        c_zx = 0.05, c_zy = 0.05,
        u0 = rand(9),
        n_maxtries = 300)

    n_tries = 0

    while n_tries <= n_maxtries
        ω₁ == nothing ? ω₁ = rand(Dω₁) : ω₁ = ω₁
        ω₂ == nothing ? ω₂ = rand(Dω₂) : ω₂ = ω₂
        ω₃ == nothing ? ω₃ = rand(Dω₃) : ω₃ = ω₃
        a₁ == nothing ? a₁ = rand(Da₁) : a₁ = a₁
        a₂ == nothing ? a₂ = rand(Da₂) : a₂ = a₂
        a₃ == nothing ? a₃ = rand(Da₃) : a₃ = a₃
        b₁ == nothing ? b₁ = rand(Db₁) : b₁ = b₁
        b₂ == nothing ? b₂ = rand(Db₂) : b₂ = b₂
        b₃ == nothing ? b₃ = rand(Db₃) : b₃ = b₃
        c₁ == nothing ? c₁ = rand(Dc₁) : c₁ = c₁
        c₂ == nothing ? c₂ = rand(Dc₂) : c₂ = c₂
        c₃ == nothing ? c₃ = rand(Dc₃) : c₃ = c₃
        pts = rossler_rossler_rossler_bidir_forced(npts,
            sample_dt = sample_dt, dt = dt, n_transient = Ttr,
            c_xy = c_xy,  c_yx = c_yx,
            c_zx = c_zx, c_zy = c_zy,
            ω₁ = ω₁, ω₂ = ω₂, ω₃ = ω₃,
            a₁ = a₁, a₂ = a₂, a₃ = a₃,
            b₁ = b₁, b₂ = b₂, b₃ = b₃,
            c₁ = c₁, c₂ = c₂, c₃ = c₃)

        if all(Matrix(pts) .< 1e10) && length(unique(pts)) > npts/2
            return pts
        end
        println("no attractor found. trying with new initial condition and parameters")
        n_tries += 1
    end
end


text:
#############################################################################

function λtoRGB(λ::T, gamma::T = 0.8) where {T<:Real}
    wavelength = λ * 1000 # λ is in um, need in nm
    if (wavelength >= 380 && wavelength <= 440)
        attenuation = 0.3 + 0.7 * (wavelength - 380) / (440 - 380)
        R = ((-(wavelength - 440) / (440 - 380)) * attenuation)^gamma
        G = 0.0
        B = (1.0 * attenuation)^gamma
    elseif (wavelength >= 440 && wavelength <= 490)
        R = 0.0
        G = ((wavelength - 440) / (490 - 440))^gamma
        B = 1.0
    elseif (wavelength >= 490 && wavelength <= 510)
        R = 0.0
        G = 1.0
        B = (-(wavelength - 510) / (510 - 490))^gamma
    elseif (wavelength >= 510 && wavelength <= 580)
        R = ((wavelength - 510) / (580 - 510))^gamma
        G = 1.0
        B = 0.0
    elseif (wavelength >= 580 && wavelength <= 645)
        R = 1.0
        G = (-(wavelength - 645) / (645 - 580))^gamma
        B = 0.0
    elseif (wavelength >= 645 && wavelength <= 750)
        attenuation = 0.3 + 0.7 * (750 - wavelength) / (750 - 645)
        R = (1.0 * attenuation)^gamma
        G = 0.0
        B = 0.0
    else
        R = 0.0
        G = 0.0
        B = 0.0
    end
    return RGB(R, G, B)
end

indexedcolor(i::Int) = ColorSchemes.hsv[rem(i / (2.1 * π), 1.0)]
indexedcolor2(i::Int) = ColorSchemes.hsv[1.0 - rem(i / (2.1 * π), 1.0)] .* 0.5

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

# main draw function
# function draw(obj; kwargs...)
#     scene, lscene = Vis.scene(resolution)
#     draw!(lscene, ob; kwargs...)
#     display(scene)

#     if (get_current_mode() == :pluto || get_current_mode() == :docs)
#         return scene
#     end
# end

# function draw!(ob; kwargs...)
#     if current_3d_scene === nothing
#         scene, lscene = Vis.scene()
#     else
#         scene = current_main_scene
#         lscene = current_3d_scene
#     end
#     draw!(lscene, ob; kwargs...)
#     display(scene)

#     if (get_current_mode() == :pluto || get_current_mode() == :docs)
#         return scene
#     end
# end

"""
    draw!(scene::Makie.LScene, csg::Union{CSGTree,CSGGenerator}; numdivisions::Int = 20, kwargs...)

Convert a CSG object ([`CSGTree`](@ref) or [`CSGGenerator`](@ref)) to a mesh using [`makemesh`](@ref) with resolution set by `numdivisions` and draw the resulting [`TriangleMesh`](@ref).
"""
draw!(scene::Scene, csg::CSGTree{T}; numdivisions::Int = 30, kwargs...) where {T<:Real} = draw!(scene, makemesh(csg, numdivisions); kwargs...)
draw!(scene::Scene, csg::CSGGenerator{T}; kwargs...) where {T<:Real} = draw!(scene, csg(); kwargs...)

#-----------------------------------------------------------------------------------------------
#   MESH from FILE
#-----------------------------------------------------------------------------------------------
function draw!(scene::Scene, ob::AbstractString; kwargs...)
    if any(endswith(lowercase(ob), x) for x in [".obj", "ply", ".2dm", ".off", ".stl"])
        meshdata = FileIO.load(ob)
        return draw!(scene, meshdata; kwargs...)
    else
        @error "Unsupported file type"
    end
end

#-----------------------------------------------------------------------------------------------
#   MESH
#-----------------------------------------------------------------------------------------------
function draw!(scene::Scene, mesh::GeometryBasics.AbstractMesh; parent_so::AbstractSceneObject = root(scene), kwargs...)
    name = "Mesh-$(UUIDs.uuid1())"
    mat = Material(;kwargs...)
    # @show typeof(mesh)
    # @show name
    scene_mesh = Mesh(mesh=mesh, material=mat, name=name)

    OpticSimVis.parent!(scene_mesh , parent_so)

    return Dict{Symbol, Any}(
        :scene => scene,
        :scene_object => scene_mesh,
    )
end

#-----------------------------------------------------------------------------------------------
#   Points
#-----------------------------------------------------------------------------------------------
function draw!(scene::Scene, points::Vector{SVector{3, T}}; parent_so::AbstractSceneObject = root(scene), kwargs...) where {T<:Real}

    name = "PointCloud-$(UUIDs.uuid1())"
    mat = Material(;kwargs...)
    pc = PointCloud(points=points, material=mat, name=name)

    OpticSimVis.parent!(pc , parent_so)

    return Dict{Symbol, Any}(
        :scene => scene,
        :scene_object => pc,
    )

    # if normals
    #     @warn "Normals being drawn from triangulated mesh, precision may be low"
    #     norigins = [Makie.Point3f0(centroid(t)) for t in tmesh.triangles[1:10:end]]
    #     ndirs = [Makie.Point3f0(normal(t)) for t in tmesh.triangles[1:10:end]]
    #     if length(norigins) > 0
    #         Makie.arrows!(scene, norigins, ndirs, arrowsize = 0.2, arrowcolor = normalcolor, linecolor = normalcolor, linewidth = 2)
    #     end
    # end
end

#-----------------------------------------------------------------------------------------------
#   TriangleMesh (OpticSim type)
#-----------------------------------------------------------------------------------------------
function draw!(scene::Scene, tmesh::TriangleMesh{T}; color=:orange, overrideAlpha = 1.0, kwargs...) where {T<:Real}
    if (overrideAlpha != 1.0)
        color = Glimmer.to_color(color)
        color = RGBA(color.r, color.g, color.b, color.alpha * overrideAlpha)
    end
    mesh = to_mesh(tmesh)
    return draw!(scene, mesh; color=color, kwargs...)
    # if normals
    #     @warn "Normals being drawn from triangulated mesh, precision may be low"
    #     norigins = [Makie.Point3f0(centroid(t)) for t in tmesh.triangles[1:10:end]]
    #     ndirs = [Makie.Point3f0(normal(t)) for t in tmesh.triangles[1:10:end]]
    #     if length(norigins) > 0
    #         Makie.arrows!(scene, norigins, ndirs, arrowsize = 0.2, arrowcolor = normalcolor, linecolor = normalcolor, linewidth = 2)
    #     end
    # end
end

#-----------------------------------------------------------------------------------------------
#   Surface{T}
#-----------------------------------------------------------------------------------------------
function draw!(scene::Scene, surf::Surface{T}; numdivisions::Int = 30, kwargs...) where {T<:Real}
    tmesh = makemesh(surf, numdivisions)
    if nothing === tmesh
        return
    end
    return draw!(scene, tmesh; kwargs...)
    # if normals
    #     ndirs = Makie.Point3f0.(samplesurface(surf, normal, numdivisions ÷ 10))
    #     norigins = Makie.Point3f0.(samplesurface(surf, point, numdivisions ÷ 10))
    #     Makie.arrows!(scene, norigins, ndirs, arrowsize = 0.2, arrowcolor = normalcolor, linecolor = normalcolor, linewidth = 2)
    # end
end

#-----------------------------------------------------------------------------------------------
# LensAssembly
#-----------------------------------------------------------------------------------------------
function draw!(scene::Scene, ass::LensAssembly{T}; drawStops = true, kwargs...) where {T<:Real}
    # @info "Draw LensAssembly"
    for (i, e) in enumerate(elements(ass))
        if (e isa OpticSim.StopSurface && drawStops || !(e isa OpticSim.StopSurface))
            draw!(scene, e; kwargs..., color = indexedcolor2(i))
        end
    end
end

#-----------------------------------------------------------------------------------------------
# sys::AbstractOpticalSystem
#-----------------------------------------------------------------------------------------------

function draw!(scene::Scene, sys::CSGOpticalSystem{T}; kwargs...) where {T<:Real}
    # @info "Draw CSGOpticalSystem"
    draw!(scene, sys.assembly; kwargs...)
    draw!(scene, sys.detector; kwargs...)
end

draw!(scene::Scene, sys::AxisymmetricOpticalSystem{T}; kwargs...) where {T<:Real} = draw!(scene, sys.system; kwargs...)

#-----------------------------------------------------------------------------------------------
#
#-----------------------------------------------------------------------------------------------
onlydetectorrays(system::Q, tracevalue::LensTrace{T,3}) where {T<:Real,Q<:AbstractOpticalSystem{T}} = onsurface(OpticSim.detector(system), point(tracevalue))

function drawtracerays!(
    scene::Scene,
    system::Q;
    raygenerator::S = Source(
        transform = translation(0.0,0.0,10.0),
        origins = Origins.RectGrid(10.0,10.0,25,25),
        directions = Constant(0.0,0.0,-1.0)
    ),
    test::Bool = false,
    trackallrays::Bool = false,
    colorbysourcenum::Bool = false,
    colorbynhits::Bool = false,
    rayfilter::Union{Nothing,Function} = onlydetectorrays,
    verbose::Bool = false,
    drawsys::Bool = false,
    drawgen::Bool = false,
    drawrays::Bool = true,
    kwargs...
) where {T<:Real,Q<:AbstractOpticalSystem{T},S<:AbstractRayGenerator{T}}

    raylines = Vector{LensTrace{T,3}}(undef, 0)

    # @info "drawtracerays!"

    if (drawgen || drawsys || drawrays)
        trace_rays_so = EmptySceneObject(name="TraceRays-$(UUIDs.uuid1())")
        OpticSimVis.parent!(trace_rays_so , root(scene))
    end

    drawgen && draw!(scene, raygenerator, norays = true, parent_so=trace_rays_so; kwargs...)
    drawsys && draw!(scene, system, parent_so=trace_rays_so; kwargs...)

    verbose && println("Tracing...")
    for (i, r) in enumerate(raygenerator)
        if i % 1000 == 0 && verbose
            print("\r $i / $(length(raygenerator))")
        end
        allrays = Vector{LensTrace{T,3}}(undef, 0)
        if trackallrays
            res = trace(system, r, trackrays = allrays, test = test)
        else
            res = trace(system, r, test = test)
        end

        if (drawrays)
            if trackallrays && !isempty(allrays)
                if rayfilter === nothing || rayfilter(system, allrays[end])
                    # filter on the trace that is hitting the detector
                    for r in allrays
                        push!(raylines, r)
                    end
                end
            elseif res !== nothing
                if rayfilter === nothing || rayfilter(system, res)
                    push!(raylines, res)
                end
            end
        end
    end
    verbose && print("\r")
    # @info "Ray Lines", length(raylines)

    if (drawrays)
        verbose && println("Drawing Rays...")
        draw!(
            scene,
            raylines,
            parent_so = trace_rays_so,
            colorbysourcenum = colorbysourcenum,
            colorbynhits = colorbynhits;
            kwargs...
        )
    end

    return raylines
end

function draw!(
    scene::Scene,
    traces::AbstractVector{LensTrace{T,N}};
    colorbysourcenum::Bool = false,
    colorbynhits::Bool = false,
    parent_so::AbstractSceneObject = root(scene),
    kwargs...
) where {T<:Real,N}
    # @info "Draw Vector of LensTrace"
    traces_by_colors = Dict()
    for trace in traces
        # draw!(scene, t; kwargs...)

        if colorbysourcenum
            color = indexedcolor(sourcenum(trace))
        elseif colorbynhits
            color = indexedcolor(OpticSim.nhits(trace))
        else
            color = λtoRGB(wavelength(trace))
        end

        final_color = RGBA(color.r, color.g, color.b, sqrt(power(trace)))

        if (!haskey(traces_by_colors, final_color))
            traces_by_colors[final_color] = [trace]
        else
            push!(traces_by_colors[final_color], trace)
        end
        # draw!(scene, (origin(ray(trace)), point(intersection(trace))); kwargs..., color = RGBA(color.r, color.g, color.b, sqrt(power(trace))), transparency = true)
    end

    # draw the lines of the rays_by_colors
    for (color, traces) in traces_by_colors
        points = Vector{SVector{3, Float64}}(undef, length(traces) * 2)
        index = 1
        for trace in traces
            p1 = origin(ray(trace))
            p2 = point(intersection(trace))
            points[index+0] = p1
            points[index+1] = p2
            index += 2
        end

        # @info typeof(points)
        segments_mat = OpticSimVis.Material(color=color)
        segments = OpticSimVis.LineSegments(
            points=points,
            material=segments_mat,
            name="Traced $(color)")
        OpticSimVis.parent!(segments, parent_so)

    end
end


text: include("basicmetrics/basicmetrics.jl")