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

using Markdown
using InteractiveUtils

# ╔═╡ ed395a10-6eef-11eb-2fa4-dfcefc661d71
using GMT

# ╔═╡ 6e509b82-6f2d-11eb-06a5-8f0c95ad5ecb
md"# Spirals
This started when Kristoff wanted to pain patterns on maps depicting sandy soils.
See the all story in this [Forum discussion](https://forum.generic-mapping-tools.org/t/how-to-generate-random-coordinates-for-pattern-fill/1322/19)
"

# ╔═╡ c95d34e0-6eef-11eb-0d6b-050d0a2904eb
md"## Archimedean spiral
From https://en.wikipedia.org/wiki/Archimedean_spiral
"

# ╔═╡ 0254aed0-6ef1-11eb-1893-55b56cbc6076
md"The begin end of this block **are NOT** the GMT's modern mode *gmt begin ... gmt end*. They are forced by the Nootebook to have more than one command per cell."

# ╔═╡ 590d7540-6ef1-11eb-2db5-316ba4efd866
begin
	# Play arround with these parameters

	T = 1;
	omega = omega = 2pi / T;
	v = 0.2;
	t = 0:0.01:5pi;
	x = v.*t .* cos.(omega .* t);
	y = v.*t .* sin.(omega .* t);
end;

# ╔═╡ 2d09b770-6ef0-11eb-3f50-0388a2cb1d46
plot(x, y, aspect=:equal, show=true)

# ╔═╡ 08699a70-6efa-11eb-1a17-1d1a287490a8
md"## Fermati spiral"

# ╔═╡ 8457d480-6efa-11eb-3a2e-f137efb6c1e8
teta = 0:0.01:5pi;

# ╔═╡ a6f62d10-6efb-11eb-3004-4dbc1a5fe7fd
xf = sqrt.(teta) .* cos.(teta);

# ╔═╡ a6f7b3b0-6efb-11eb-2004-8b152d668695
yf = sqrt.(teta) .* sin.(teta);

# ╔═╡ a6f87700-6efb-11eb-2443-3b78c63192b9
plot(xf,yf, aspect=:equal, show=true)

# ╔═╡ 04a8cc70-6efb-11eb-3aea-df71f7e7b780
md"## Sunflower
From [This FEX contribution](https://www.mathworks.com/matlabcentral/fileexchange/10796-model-a-sunflower-with-the-golden-ratio). The author here wanted to reflect the fact that on a sunflower the seeds close to the center are smaller and have a higher density.
"

# ╔═╡ 1842f4e0-6efb-11eb-17dd-c74c91e5792f
phi = (sqrt(5)-1)/2;

# ╔═╡ 1c958620-6efb-11eb-35f9-d51c907524fe
n = 2618;

# ╔═╡ 38c2cfb0-6efb-11eb-1daf-c1883c2b8443
rho = (2:n-1) .^ phi;

# ╔═╡ 802977a0-6efb-11eb-1f34-997996793d9e
theta = (2:n-1)*2pi*phi;

# ╔═╡ 8c5ca2e0-6efb-11eb-2e5b-3fe95be4a166
scatter(rho .* cos.(theta), rho .* sin.(theta), marker=:point, aspect=:equal, show=true)

# ╔═╡ d5e2e7f0-6f26-11eb-3bc0-2bb0532670ff
md"## Another Sunflower
This one was reversed from the javascript in [this page](https://github.com/jacquerie/sunflower/blob/gh-pages/js/application.js), which follows the original work of Helmut Vogel in [A better Way to Construct the Sunflower Head](http://dx.doi.org/10.1016%2F0025-5564%2879%2990080-4), where he proposed that spiral branches of seeds in a sunflower head are added from the center at an angle of 137.5∘ from the preceding one.

This time we will also color the seed points in function of *r*, the distance to the center and pain with a dark background.
"

# ╔═╡ 3ed0c890-6f27-11eb-0395-d309ad06bd06
begin
	angle = 137.5;	# Play with this angle between [137.0 138.0]. Amazing the effect, no?
	alfa = 2pi * angle / 360;
	n_seeds = 1500;
	seeds = 0:n_seeds;
	r = sqrt.(seeds);
	ϕ = alfa * seeds;
	C = makecpt(range=(1,sqrt(n_seeds),1), cmap=:buda);	# Color map to paint the seeds
	scatter(r .* cos.(ϕ), r .* sin.(ϕ), marker=:point, cmap=C, zcolor=r,
	#scatter(r, ϕ, proj=:polar, marker=:point, cmap=C, zcolor=r,
		frame=(fill=20,), aspect=:equal, show=true)
end

# ╔═╡ Cell order:
# ╟─6e509b82-6f2d-11eb-06a5-8f0c95ad5ecb
# ╟─c95d34e0-6eef-11eb-0d6b-050d0a2904eb
# ╠═ed395a10-6eef-11eb-2fa4-dfcefc661d71
# ╟─0254aed0-6ef1-11eb-1893-55b56cbc6076
# ╠═590d7540-6ef1-11eb-2db5-316ba4efd866
# ╠═2d09b770-6ef0-11eb-3f50-0388a2cb1d46
# ╠═08699a70-6efa-11eb-1a17-1d1a287490a8
# ╠═8457d480-6efa-11eb-3a2e-f137efb6c1e8
# ╠═a6f62d10-6efb-11eb-3004-4dbc1a5fe7fd
# ╠═a6f7b3b0-6efb-11eb-2004-8b152d668695
# ╠═a6f87700-6efb-11eb-2443-3b78c63192b9
# ╟─04a8cc70-6efb-11eb-3aea-df71f7e7b780
# ╠═1842f4e0-6efb-11eb-17dd-c74c91e5792f
# ╠═1c958620-6efb-11eb-35f9-d51c907524fe
# ╠═38c2cfb0-6efb-11eb-1daf-c1883c2b8443
# ╠═802977a0-6efb-11eb-1f34-997996793d9e
# ╠═8c5ca2e0-6efb-11eb-2e5b-3fe95be4a166
# ╟─d5e2e7f0-6f26-11eb-3bc0-2bb0532670ff
# ╠═3ed0c890-6f27-11eb-0395-d309ad06bd06


text:
# julia Cut.jl --parallel no --lmax 5
##################################################

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

aMeasure = 10.0 # Semi-major axis (in mpc) along which the coefficients Djj are computed

##################################################
using HDF5 # To have access to .hf5 files
########################################
include("../sources/Main.jl") # Loading the main code
########################################
jminMeasure, jmaxMeasure = 0.01,0.999 # Range in j where the Djj are computed
nbjMeasure = 100 # Number of j for which the Djj are computed
tabjMeasure = exp.(range(log(jminMeasure),length=nbjMeasure,log(jmaxMeasure)))

const tabDRRjj = zeros(Float64,nbjMeasure) # Values of the DRR_jj coefficients on the (a,j)-grid
const tabDNRjj = zeros(Float64,nbjMeasure) # Values of the DRR_jj coefficients on the (a,j)-grid

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

function tabDRRjj!()
    if (PARALLEL == "yes") # Computation is made in parallel
        Threads.@threads for ij=1:nbjMeasure # Loop over the angular momenta
            jMeasure = tabjMeasure[ij] # Current angular momentum
            if jMeasure<jlc(aMeasure) # Do not compute DRR within the loss cone
                tabDRRjj[ij] = 0.0
            else
                tabSMARes_parallel = zeros(Float64,2,2) # Container of tabSMARes for the current thread
                IntTable_K_parallel = IntTable_create!()

                DRR = DRR_jj(aMeasure,jMeasure,IntTable_K_parallel,tabSMARes_parallel)
                tabDRRjj[ij] = DRR # Computing DRR_jj with the thread's containers
            end
        end
    else # Computation is not made in parallel
        for ij=1:nbjMeasure # Loop over the angular momenta
            jMeasure = tabjMeasure[ij] # Current angular momentum
            if jMeasure<jlc(aMeasure) # Do not compute DRR within the loss cone
                tabDRRjj[ij] = 0.0
            else
                DRR = DRR_jj(aMeasure,jMeasure)
                tabDRRjj[ij] = DRR # Computing DRR_jj with the global containers
            end
        end
    end
end

function tabDNRjj!()
    if (PARALLEL == "yes") # Computation is made in parallel
        Threads.@threads for ij=1:nbjMeasure # Loop over the angular momenta
            jMeasure = tabjMeasure[ij] # Current angular momentum
            if jMeasure<jlc(aMeasure) # Do not compute DNR within the loss cone
                tabDNRjj[ij] = 0.0
            else
                DNR = DNR_jj(aMeasure,jMeasure)
                tabDNRjj[ij] = DNR # Computing DNR_jj with the thread's containers
            end
        end
    else # Computation is not made in parallel
        for ij=1:nbjMeasure # Loop over the angular momenta
            jMeasure = tabjMeasure[ij] # Current angular momentum
            if jMeasure<jlc(aMeasure) # Do not compute DNR within the loss cone
                tabDNRjj[ij] = 0.0
            else
                DNR = DNR_jj(aMeasure,jMeasure)
                tabDNRjj[ij] = DNR # Computing DNR_jj with the global containers
            end
        end
    end
end

########################################
namefile = "../data/Dump_Diffusion_Coefficients_Cut.hf5"
########################################
# Function that writes to .hf5 files
########################################
function writedump!(namefile)
    file = h5open(namefile,"w") # Opening the file
    write(file,"tabj",tabjMeasure) # Dumping the list of angular momenta
    write(file,"tabDRRjj",tabDRRjj) # Dumping the DRRjj(a,j) for the current used value of lmax
    write(file,"tabDNRjj",tabDNRjj) # Dumping the DNRjj(a,j) for the current used value of lmax
    write(file,"aMeasure",aMeasure) # Dumping aMeasure
    write(file,"jmin",jminMeasure) # Dumping jminMeasure
    write(file,"jmax",jmaxMeasure) # Dumping jmaxMeasure
    write(file,"nbj",nbjMeasure) # Dumping nbjMeasure
    write(file,"jlc",jlc(aMeasure)) # Dumping jlc(a)
    close(file) # Closing the file
end

########################################
println("Computing SRR diffusion coefficients...")
@time tabDRRjj!()
println("Computing NR diffusion coefficients...")
@time tabDNRjj!()

########################################
println("Saving data...")
writedump!(namefile) # Dumping the computed table

ext: jl


text:
export VRepOptimizer

"""
    VRepOptimizer{T} <: AbstractPolyhedraOptimizer{T}

Linear Programming solver using the V-representation of the feasible set to find
the optimal solution.
"""
mutable struct VRepOptimizer{T} <: AbstractPolyhedraOptimizer{T}
    library::Union{Nothing, Library}
    lphrep::LPHRep{T}
    # Feasible set defined by user
    rep::Union{Rep{T}, Nothing}
    # Either `rep` or `polyhedron(lphrep, library)`.
    # It is kept between consecutive solve if not modified,
    # e.g. if only the objective is changed.
    feasible_set::Union{Rep{T}, Nothing}

    objective_sense::MOI.OptimizationSense
    objective_func::Union{SparseVector{T, Int64}, Nothing}
    objective_constant::T

    status::MOI.TerminationStatusCode
    solution::Union{AbstractVector{T}, Nothing}

    function VRepOptimizer{T}(library::Union{Nothing, Library} = nothing) where T
        new(library, LPHRep(_MOIModel{T}()), nothing, nothing,
            MOI.FEASIBILITY_SENSE, nothing, zero(T),
            MOI.OPTIMIZE_NOT_CALLED, nothing)
    end
end

coefficient_type(::VRepOptimizer{T}) where {T} = T
MOI.get(::VRepOptimizer, ::MOI.SolverName) = "VRep"

function MOI.empty!(lpm::VRepOptimizer{T}) where T
    lpm.lphrep = LPHRep(_MOIModel{T}())
    lpm.rep = nothing
    lpm.feasible_set = nothing
    lpm.objective_sense = MOI.FEASIBILITY_SENSE
    lpm.objective_func = nothing
    lpm.objective_constant = zero(T)
    lpm.status = MOI.OPTIMIZE_NOT_CALLED
    lpm.solution = nothing
end
function MOI.is_empty(lpm::VRepOptimizer{T}) where T
    MOI.is_empty(lpm.lphrep.model) &&
    lpm.rep === nothing &&
    lpm.feasible_set === nothing &&
    lpm.objective_sense == MOI.FEASIBILITY_SENSE &&
    lpm.objective_func === nothing &&
    iszero(lpm.objective_constant) &&
    lpm.status == MOI.OPTIMIZE_NOT_CALLED &&
    lpm.solution === nothing
end

function MOI.optimize!(lpm::VRepOptimizer{T}) where T
    if lpm.rep === nothing
        lpm.feasible_set = lpm.lphrep
    else
        if hasallhalfspaces(lpm.lphrep)
            error("Cannot provide both a polyhedral feasible set and additional constraints.")
        end
        lpm.feasible_set = lpm.rep
    end
    if lpm.feasible_set isa HRepresentation
        if lpm.library === nothing
            lpm.feasible_set = polyhedron(lpm.feasible_set)
        else
            lpm.feasible_set = polyhedron(lpm.feasible_set, lpm.library)
        end
    end
    if lpm.feasible_set isa VRepresentation
        prob = lpm.feasible_set
    else
        @assert lpm.feasible_set isa Polyhedron
        prob = vrep(lpm.feasible_set)
    end
    N = fulldim(prob)
    if !haspoints(prob) && !haslines(prob) && !hasrays(prob)
        lpm.status = MOI.INFEASIBLE
        lpm.solution = nothing
    elseif lpm.objective_sense == MOI.FEASIBILITY_SENSE
        lpm.status = MOI.OPTIMAL
        lpm.solution = first(points(prob))
    else
        better(a, b) = (lpm.objective_sense == MOI.MAX_SENSE ? a > b : a < b)
        _better(a, b) = (lpm.objective_sense == MOI.MAX_SENSE ? _gt(a, b) : _lt(a, b))
        bestobjval = zero(T)
        lpm.solution = nothing
        for r in allrays(prob)
            objval = lpm.objective_func ⋅ r
            if _better(objval, bestobjval)
                bestobjval = objval
                lpm.solution = coord(r)
            end
        end
        if lpm.solution !== nothing
            lpm.status = MOI.DUAL_INFEASIBLE
        else
            for p in points(prob)
                objval = lpm.objective_func ⋅ p
                if lpm.solution === nothing || better(objval, bestobjval)
                    bestobjval = objval
                    lpm.solution = p
                end
            end
            lpm.status = MOI.OPTIMAL
        end
        @assert lpm.solution !== nothing
    end
end

MOI.get(lpm::VRepOptimizer, ::MOI.TerminationStatus) = lpm.status
function MOI.get(lpm::VRepOptimizer, ::MOI.ResultCount)
    if lpm.status == MOI.OPTIMAL || lpm.status == MOI.DUAL_INFEASIBLE
        return 1
    else
        return 0
    end
end
function MOI.get(lpm::VRepOptimizer{T}, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},
                                         <:Union{MOI.EqualTo{T},
                                                 MOI.LessThan{T}}}) where T
    return MOIU.get_fallback(lpm, attr, ci)
end
function MOI.get(lpm::VRepOptimizer, ::MOI.ObjectiveValue)
    if lpm.status == MOI.OPTIMAL
        return lpm.objective_func ⋅ lpm.solution + lpm.objective_constant
    elseif lpm.status == MOI.DUAL_INFEASIBLE
        return lpm.objective_func ⋅ lpm.solution
    else
        error("No objective value available when termination status is $(lpm.status).")
    end
end
function MOI.get(lpm::VRepOptimizer, ::MOI.PrimalStatus)
    if lpm.status == MOI.OPTIMAL
        return MOI.FEASIBLE_POINT
    elseif lpm.status == MOI.DUAL_INFEASIBLE
        return MOI.INFEASIBILITY_CERTIFICATE
    else
        return MOI.NO_SOLUTION
    end
end
function MOI.get(lpm::VRepOptimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex)
    if lpm.status != MOI.OPTIMAL && lpm.status != MOI.DUAL_INFEASIBLE
        error("No primal value available when termination status is $(lpm.status).")
    end
    return lpm.solution[vi.value]
end


text:
#
# see: https://github.com/JuliaWeb/Mux.jl
# and:
#
using Mux
import Mux.WebSockets
using JSON
using HttpCommon
using Logging, LoggingExtras

using BudgetConstraints
using PovertyAndInequalityMeasures

import ScottishTaxBenefitModel: Definitions, FRSHouseholdGetter, ExampleHouseholdGetter
   ModelHousehold, Utils, MiniTB, GeneralTaxComponents, WebModelLibs
using .Definitions
using .FRSHouseholdGetter
using .ExampleHouseholdGetter
using .ModelHousehold
using .Utils
using .MiniTB
using .GeneralTaxComponents
using .WebModelLibs

## !!! needs Julia 1.3:w
import Base.Threads.@spawn

const DEFAULT_PORT=8000
const DEFAULT_SERVER="http://localhost:$DEFAULT_PORT/"

const DEFAULT_TEST_URL="$(DEFAULT_SERVER)/bc?it_allow=300.0&it_rate_1=0.25&it_rate_2=0.5&it_band=10000&benefit1=150.0&benefit2=60.0&ben2_min_hours = 150.0&ben2_taper=0.5&ben2_u_limit = 250.0"
const ZERO_TEST_URL="$(DEFAULT_SERVER)/bc?it_allow=0&it_rate_1=0&it_rate_2=0&it_band=0&benefit1=0&benefit2=0.0&ben2_taper=0&ben2_min_hours=0&ben2_u_limit=0"

@debug "server starting up"

# gks work round some (version?) thing in the diffeq code
function get_thing( thing::AbstractArray, key :: AbstractString, default :: AbstractString )
   for i in thing
      if i[1] == key
         return i[2]
      end
   end # loop
   default
end #get

function web_get_hh( hdstr::AbstractString )
   hid = parse(Int64, hdstr)
   JSON.json(FRSHouseholdGetter.get_household(hid))
end

function addqstrdict(app, req  :: Dict )
   req[:parsed_querystring] = qstrtodict(req[:query])
   return app(req)
end

# Better error handling
function errorCatch( app, req  :: Dict )
   try
      app(req)
   catch e
      @error "Error occured!"
      io = IOBuffer()
      showerror(io, e)
      err_text = takebuf_string(io)
      @error err_text
      resp = withHeaders(JSON.json(Dict("message" => err_text, "error" => true)), req)
      resp[:status] = 500
      return resp
   end
end

function d100( v :: Number ) :: Number
   v/100.0
end

function web_map_params( req  :: Dict, defaults = MiniTB.DEFAULT_PARAMS ) :: MiniTB.TBParameters
   querydict = req[:parsed_querystring]
   tbparams = deepcopy( defaults )
   tbparams.it_allow = get_if_set("it_allow", querydict, tbparams.it_allow, operation=weeklyise )
   tbparams.it_rate[1] = get_if_set("it_rate_1", querydict, tbparams.it_rate[1], operation=d100 )
   tbparams.it_rate[2] = get_if_set("it_rate_2", querydict, tbparams.it_rate[2], operation=d100 )
   tbparams.it_band[1] = get_if_set("it_band", querydict, tbparams.it_band[1], operation=weeklyise)
   tbparams.benefit1 = get_if_set("benefit1", querydict, tbparams.benefit1)
   tbparams.benefit2 = get_if_set("benefit2", querydict, tbparams.benefit2)
   tbparams.ben2_min_hours = get_if_set("ben2_min_hours", querydict, tbparams.ben2_min_hours)
   tbparams.ben2_taper = get_if_set("ben2_taper", querydict, tbparams.ben2_taper, operation=d100)
   tbparams.ben2_u_limit = get_if_set("ben2_u_limit", querydict, tbparams.ben2_u_limit)
   tbparams.basic_income = get_if_set("basic_income", querydict, tbparams.basic_income)
   @debug "DEFAULT_PARAMS\n$DEFAULT_PARAMS"
   @debug "tbparams\n$tbparams"
   tbparams
end

const example_names, num_households, num_people = load_data( load_examples = true, load_main = true, start_year = 2015 )

const DEFAULT_BC = local_makebc(MiniTB.DEFAULT_PERSON, MiniTB.DEFAULT_PARAMS)

const BASE_RESULTS = create_base_results( num_households, num_people )

function makeztparams() :: MiniTB.TBParameters
   pars = deepcopy( MiniTB.DEFAULT_PARAMS )
   pars.it_rate = [0.0,0.0]
   pars
end

const BC_500_SETTINGS = BCSettings( maxgross=500 )

const ZERO_TAX_PARAMS = makeztparams()
const ZERO_DEFAULT_BC = local_makebc(MiniTB.DEFAULT_PERSON, MiniTB.ZERO_PARAMS)
const ZERO_TAX_BC = local_makebc(MiniTB.DEFAULT_PERSON, ZERO_TAX_PARAMS, BC_500_SETTINGS )

function main_run_to_json( tbparams :: MiniTB.TBParameters ):: String
   results = do_one_run( tbparams, num_households, num_people, NUM_REPEATS )
   summary_output = summarise_results!( results=results, base_results=BASE_RESULTS )
   JSON.json( summary_output )
end

function web_do_one_run( req :: Dict ) :: AbstractString
   @info "web_do_one_run; running on thread $(Threads.threadid())"
   tbparams = web_map_params( req )
   json = main_run_to_json( tbparams )
   # headers could include (e.g.) a timestamp, so add after caching
end # do_one_run

function web_makebc( req  :: Dict ) :: AbstractString
   @debug "web_makebc; running on thread $(Threads.threadid())"
   tbparams = web_map_params( req )
   bc =  local_makebc( DEFAULT_PERSON, tbparams )
   JSON.json((base = DEFAULT_BC, changed = bc))
end

function web_makezbc( req  :: Dict ) :: AbstractString
   @debug "web_makezbc; running on thread $(Threads.threadid())"
   tbparams = web_map_params( req, MiniTB.ZERO_PARAMS )
   bc =  local_makebc( DEFAULT_PERSON, tbparams )
   JSON.json((base = ZERO_DEFAULT_BC, changed = bc))
end

function web_makeztbc( req  :: Dict ) :: AbstractString
   @debug "web_makeztbc; running on thread $(Threads.threadid())"
   tbparams = web_map_params( req, ZERO_TAX_PARAMS )
   bc =  local_makebc( DEFAULT_PERSON, tbparams, BC_500_SETTINGS )
   JSON.json((base = ZERO_TAX_BC, changed = bc))
end

function web_doineq( req  :: Dict ) :: AbstractString
   @debug "web_doineq; running on thread $(Threads.threadid())"
   querydict = req[:parsed_querystring]
   inc=zeros(0)
   pop=zeros(0)
   for i in 1:10 # can's see how to pass a json array from js to this, so ...
      push!( inc, Float64(get_if_set("inc_$i", querydict, 0 )))
      push!( pop, Float64(get_if_set("pop_$i", querydict, 0 )))
   end
   data = hcat( pop, inc )
   ineq = PovertyAndInequalityMeasures.make_inequality( data, 1, 2 )
   JSON.json( ( data=data, ineq=ineq ))
end

# Headers -- set Access-Control-Allow-Origin for either dev or prod
# this is from https://github.com/JuliaDiffEq/DiffEqOnlineServer
#
function add_headers( json :: AbstractString ) :: Dict
    headers  = HttpCommon.headers()
    headers["Content-Type"] = "application/json; charset=utf-8"
    headers["Access-Control-Allow-Origin"] = "*"
    Dict(
       :headers => headers,
       :body=> json
    )
end

#
# This is my attempt at starting a task using the 1.3 @spawn macro
# 1st parameter is a function that returns String (probably a json string) and accepts the req Dict
#
function do_in_thread( the_func::Function, req :: Dict ) :: Dict
   response = @spawn the_func( req )
   # note that the func returns a string but response is a Future type
   # line below converts response to a string
   @debug "do_in_thread response is $response"
   json = fetch( response )
   add_headers( json )
end

# configure logger; see: https://docs.julialang.org/en/v1/stdlib/Logging/index.html
# and: https://github.com/oxinabox/LoggingExtras.jl
logger = FileLogger("/var/tmp/stb_log.txt")
global_logger(logger)
LogLevel( Logging.Info )

#
# from diffeq thingy instead of Mux.defaults
#
# ourstack = stack(Mux.todict, errorCatch, Mux.splitquery, Mux.toresponse) # from DiffEq
#
@app stbdemo = (
   Mux.defaults,
   addqstrdict,
   page( respond("
A Scottish Tax-Benefit Model

; very early version
")),
   page( "/hhld/:hid", req -> web_get_hh((req[:params][:hid]))), # note no headers
   page("/ineq", req -> do_in_thread( web_doineq, req )),
   page("/bc", req -> do_in_thread( web_makebc, req )),
   page("/zbc", req -> do_in_thread( web_makezbc, req )),
   page("/ztbc", req -> do_in_thread( web_makeztbc, req )),
   # handle main stb run a bit differently so we can
   # check & cache the results.
   # page("/stb", req -> web_do_one_run_cached( req )),
   page("/stb", req -> do_in_thread( web_do_one_run, req )),
   Mux.notfound(),
)

port = DEFAULT_PORT
if length(ARGS) > 0
   port = parse(Int64, ARGS[1])
end

serve(stbdemo, port)

while true # FIXME better way?
   @debug "main loop; server running on port $port"
   sleep( 60 )
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function _find_neworigin(p::AbstractNode, q::AbstractNode, oldorigin::AbstractNode, neworigin::AbstractNode)
	if n_neighbour(q) > 2
		neworigin = q
	else
		for link in q.links
			link.to == p && continue
			neworigin = _find_neworigin(q, link.to, oldorigin, neworigin)
			neworigin ≠ oldorigin && break
		end
	end

	return neworigin
end

"""
	pluck_nonsplitting!(tree)

Plucks out all the non-splitting nodes in a `tree`.

If the current origin of the tree is a non-splitting node, the next splitting node (in preorder) will be designated as the new origin of the tree.
"""
function pluck_nonsplitting!(tree::AbstractTree)
	targets = find_nonsplitting(tree)
	isempty(targets) && return nothing

	neworigin = tree.origin
	if targets[1] == tree.origin
		neworigin = _find_neworigin(tree.origin, tree.origin, tree.origin, tree.origin)
	end

	for p in targets
		q, r = neighbours(p)
		pluck!(p, q, r)
	end

	tree.origin = neworigin

	return nothing
end


text:
"""
read_seq_file(filename, barcode, stagger, read_length; verbose=false)

"""
function read_seq_file(mapping::Dict{DNASequence, String},
                       filename::String,
                       barcode::DNASequence,
                       stagger::Int,
                       read_length::Int,
                       truncate::Int;
                       verbose=false)
    # inflate streamed data from zipped file
    stream = open(filename) |> ZlibInflateInputStream
    try
        reader = FASTQReader(stream, quality_encoding=:illumina18)
        record = FASTQSeqRecord{DNASequence}()

        searchseq = barcode * dna"TCTTGTGGAAAGGACGAAACACCG"
        len_search = length(searchseq)
        # for detecting alternative staggers
        search_query = ExactSearchQuery(searchseq)

        sample = ScreenSample(filename, read_length, 41)

        while !eof(reader)
            read!(reader, record)
            sample.reads += 1
            # number of errors in the constant region
            dist = hamming(searchseq, record.seq[(stagger+1):(stagger+len_search)])
            # get quality scores
            qual = record.metadata.quality

            priming_range = search(record.seq, search_query)
            (priming_range.stop > 0) && (sample.alt_staggers[priming_range.start] += 1)

            if verbose && sample.reads < 100
                println("\n\n", record.seq, " Read")
                println("*"^stagger, searchseq, "*"^(read_length - len_search - stagger), " Target constant region")
                println("-"^stagger, record.seq[stagger+1:stagger+len_search], "-"^(read_length - len_search - stagger), " Read constant region")
                println("-"^(stagger+len_search), record.seq[stagger+len_search+1:end], " Read sgrna")
                println(length(record.seq[stagger+len_search+1:end]))
            end

            if dist <= 5
                # read must have a PHRED quality above 37 for entire barcode region
                sgrna_quality = qual[(stagger + len_search+1):(stagger+len_search + truncate)]

                if sum(sgrna_quality .< 37) == 0
                    for pos in 1:read_length
                        sample.hit_qualities[qual[pos], pos] += 1
                    end

                    sgrna = record.seq[(stagger + len_search + 1):(stagger + len_search + truncate)]

                    if haskey(mapping, sgrna)
                        sample.hits += 1
                        (verbose && sample.reads < 100) && println("\nMatched as $(mapping[sgrna])")
                        if haskey(sample.counts, mapping[sgrna])
                            sample.counts[mapping[sgrna]] += 1
                        else
                            sample.counts[mapping[sgrna]] = 1
                        end

                    else
                        (verbose && sample.reads < 100) && println("\nNo match for $sgrna")
                        sample.no_match += 1
                    end
                elseif verbose && sample.reads < 100
                    println("\nDoes not pass quality threshold in sgrna\n$(join(sgrna_quality, ','))")
                else
                    sample.low_quality += 1
                end
            elseif verbose && sample.reads < 100
                println("\nToo many errors in constant region ($(dist) errors detected)")
            else
                sample.constant_error += 1
            end

            for pos in 1:read_length
                sample.all_qualities[qual[pos], pos] += 1
            end
        end
        # compare the provided stagger versus the most common detected one and insure that they
        # are the same. this is a sanity check that we're not throwing away too many reads.
        alt_stagger = findmax(sample.alt_staggers)[2] - 1
        (alt_stagger != stagger) && warn("An alternative stagger of $alt_stagger was detected as the dominant stagger versus the provided stagger of $stagger.")
        return sample
    finally
        close(stream)
        end
end


text:
function randomSphereFilling(volume,r,seed)

    #Create emptyobjects
    frontSpheresQueue = []
    assemblyList = []
    gridIndex = []
    prevRejectedQueue = []
    newlyRejectedQueue = []
    δbox = 3*r

    #Create seed spheres
    ##First
    push!(frontSpheresQueue,seed)
    push!(assemblyList,seed)
    push!(gridIndex,seed)
    ##Second
    v = rand(Normal(0,1),3)
    v ./= sqrt(sum(v.^2))
    pos2 = seed .+ 2*r .*v
    push!(frontSpheresQueue,pos2)
    push!(assemblyList,pos2)
    push!(gridIndex,pos2)
    ##Third
    mid = seed .+ r .*v
    if v[1] != 0 || v[2] != 0 #Make first orthogonal vector
        orth1 = [v[2],-v[1],0]
    else
        orth1 = [0,v[3],-v[2]]
    end
    orth1 ./= sqrt(sum(orth1.^2))
    orth2 = cross(v,orth1)
    v = rand(Normal(0,1)).*orth1 + rand(Normal(0,1)).*orth2
    v ./= sqrt(sum(v.^2))
    pos2 = mid .+ 2*r*cos(pi/6) .*v
    push!(frontSpheresQueue,pos2)
    push!(assemblyList,pos2)
    push!(gridIndex,pos2)

    #while loop
    while !isempty(frontSpheresQueue)
        #Retrieve one active sphere
        currentSphere = frontSpheresQueue[1]
        #Add neighbour Spheres to currect sphere
        neighboringSpheres = []
        for s in gridIndex
            if sqrt(sum((s-currentSphere).^2)) <= δbox
                push!(neighboringSpheres,s)
            end
        end
        #Candidate point list

    end

    return
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
using brainflow

# specify logging library to use
brainflow.enable_dev_brainflow_logger(brainflow.BOARD_CONTROLLER)

params = BrainFlowInputParams()
board_shim = brainflow.BoardShim(brainflow.SYNTHETIC_BOARD, params)
sampling_rate = brainflow.get_sampling_rate(brainflow.SYNTHETIC_BOARD)
nfft = brainflow.get_nearest_power_of_two(sampling_rate)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_board_data(board_shim)
brainflow.release_session(board_shim)

eeg_channels = brainflow.get_eeg_channels(brainflow.SYNTHETIC_BOARD)
# second channel of synthetic board is sine wave at 10 Hz, should see huge 'alpha'
data_second_channel = data[eeg_channels[2], :]

# optional: detrend
brainflow.detrend(data_second_channel, brainflow.LINEAR)
# psd is a tuple of ampls and freqs
psd = brainflow.get_psd_welch(data_second_channel, nfft, Integer(nfft / 2), sampling_rate, brainflow.BLACKMAN_HARRIS)
band_power_alpha = brainflow.get_band_power(psd, 7.0, 13.0)
band_power_beta = brainflow.get_band_power(psd, 14.0, 30.0)
println(band_power_alpha / band_power_beta)

max_stars_repo_name: matthijscox/brainflow

max_issues_repo_name: matthijscox/brainflow

max_forks_repo_name: matthijscox/brainflow


text:
#Distributions.Normal()
#Distributions.TruncatedNormal()
#rand()
#Distributions.logpdf(d, x)
#Distributions.pdf(d, x)

#sample_kernel()
#density_kernel()

"""
Simultaneously generate sampler and density function for kernel
"""

function make_normal_kernel()
     function sample_kernel(x0::Float64, kernelbw::Float64)
        return randn() * kernelbw + x0
    end
     function density_kernel(x::Float64, x0::Float64, kernelbw::Float64)
        d = Distributions.Normal(x0, kernelbw)
        return Distributions.pdf(d, x)
     end
    return sample_kernel, density_kernel
end

function make_jointnormal_kernel()
     function sample_kernel(x0::Array{Float64, 1}, kernelbw::Array{Float64,1})
        return randn(length(x0)) .* kernelbw + x0
    end
    function density_kernel(x::Array{Float64, 1}, x0::Array{Float64, 1}, kernelbw::Array{Float64, 1})
        dA = [Distributions.Normal(x0i, ksd) for (x0i, ksd) in zip(x0, kernelbw)]
        logpA = [Distributions.logpdf(d, xi) for (d, xi) in zip(dA, x)]
        logp = sum(logpA)
        return exp(logp)
    end
    return sample_kernel, density_kernel
end

function make_truncatednormal_kernel(lowerbound::Float64, upperbound::Float64)
    function sample_kernel(x0::Float64, kernelbw::Float64)
        d = Distributions.TruncatedNormal(x0, kernelbw, lowerbound, upperbound)
        return rand(d)
    end
    function density_kernel(x::Float64, x0::Float64, kernelbw::Float64)
        d = Distributions.TruncatedNormal(x0, kernelbw, lowerbound, upperbound)
        return Distributions.pdf(d, x)
    end
    return sample_kernel, density_kernel
end

function make_truncatednormal_kernel(lowerbounds::Array{Float64, 1}, upperbounds::Array{Float64, 1})
    function sample_kernel(x0::Array{Float64, 1}, kernelbw::Array{Float64, 1})
        dA = [Distributions.TruncatedNormal(x0i, ksd, l, u) for (x0i, ksd, l, u) in zip(x0, kernelbw, lowerbounds, upperbounds)]
        return rand.(dA)
    end
    function density_kernel(x::Array{Float64, 1}, x0::Array{Float64, 1}, kernelbw::Array{Float64, 1})
        dA = [Distributions.TruncatedNormal(x0i, ksd, l, u) for (x0i, ksd, l, u) in zip(x0, kernelbw, lowerbounds, upperbounds)]
        logpA = [Distributions.logpdf(d, xi) for (d, xi) in zip(dA, x)]
        logp = sum(logpA)
        return exp(logp)
    end
    return sample_kernel, density_kernel
end

function make_truncatednormal_kernel(bounds::Array{Tuple{Float64, Float64}, 1})
    function sample_kernel(x0::Array{Float64, 1}, kernelbw::Array{Float64, 1})
        dA = [Distributions.TruncatedNormal(x0i, ksd, b...) for (x0i, ksd, b) in zip(x0, kernelbw, bounds)]
        return rand.(dA)
    end
    function density_kernel(x::Array{Float64, 1}, x0::Array{Float64, 1}, kernelbw::Array{Float64, 1})
        dA = [Distributions.TruncatedNormal(x0i, ksd, b...) for (x0i, ksd, b) in zip(x0, kernelbw, bounds)]
        logpA = [Distributions.logpdf(d, xi) for (d, xi) in zip(dA, x)]
        logp = sum(logpA)
        return exp(logp)
    end
    return sample_kernel, density_kernel
end

function make_joint_kernel(distr::Array{Any, 1}, params::Array{Tuple, 1})
    function sample_kernel(x0::Array{Float64, 1}, kernelbw::Array{Float64, 1})
        dA = [d(x0i, ksd, p...) for (d, x0i, ksd, p) in zip(distr, x0, kernelbw, params)]
        return rand.(dA)
    end
    function density_kernel(x::Array{Float64, 1}, x0::Array{Float64, 1}, kernelbw::Array{Float64, 1})
        dA = [d(x0i, ksd, p...) for (d, x0i, ksd, p) in zip(distr, x0, kernelbw, params)]
        logpA = [Distributions.logpdf(d, xi) for (d, xi) in zip(dA, x)]
        logp = sum(logpA)
        return exp(logp)
    end
    return sample_kernel, density_kernel
end


text:
__precompile__(true)

module EarthEngine

using Logging
using PyCall
using OrderedCollections
import Base:
    length,
    keys,
    contains,
    split,
    replace,
    lowercase,
    string,
    filter,
    union,
    size,
    identity,
    map,
    first,
    get,
    repeat

const module_dir = @__DIR__

# define a constant pointer to the ee object
const ee = PyNULL()

version() = VersionNumber(ee.__version__)

const pre_type_map = []

# Maps a python object corresponding to an EE class to a Julia type which
# wraps that class.
const type_map = OrderedDict{PyObject,DataType}()

abstract type EEWrapped end

PyCall.PyObject(x::EEWrapped) = x.pyo

"""
    @pytype name class

Macro for creating a Julia Type that wraps a PyObject class
"""
macro pytype(name, class)
    quote
        struct $(name) <: EEWrapped
            pyo::PyObject
            $(esc(name))(pyo::PyObject) = new(pyo)
            function $(esc(name))(args...; kwargs...)
                ee_method = ($class)()
                new(pycall(ee_method, PyObject, args...; kwargs...))
            end
        end

        # This won't work until PyCall is updated to support
        # the Julia 1.0 iteration protocol.
        function Base.iterate(x::$name, state...)
            res = Base.iterate(x.pyo, state...)
            if res === nothing
                return nothing
            else
                value, state = res
                return ee_wrap(value), state
            end
        end

        push!(pre_type_map, ($class, $name))
    end
end

"""
    ee_wrap(pyo::PyObject)

Function for wrapping a Python object defined in the type map
"""
function ee_wrap(pyo::PyObject)
    for (pyt, pyv) in type_map
        pyt === nothing && continue
        if pyisinstance(pyo, pyt)
            return pyv(pyo)
        end
    end
    return convert(PyAny, pyo)
end

quot(x) = Expr(:quote, x)

ee_wrap(x::Union{AbstractArray,Tuple}) = [ee_wrap(_) for _ in x]

ee_wrap(pyo) = pyo

fix_arg(x) = x

pyattr(class, method) = pyattr(class, method, method)

"""
    pyattr(class method orig_method)

Function for creating a method signature for a Julia Type
In Python world method(T::class) is analagous to class.method()
"""
function pyattr(class, jl_method, py_method)
    quote
        function $(esc(jl_method))(pyt::$class, args...; kwargs...)
            new_args = fix_arg.(args)
            method = pyt.pyo.$(string(py_method))
            pyo = pycall(method, PyObject, new_args...; kwargs...)
            wrapped = ee_wrap(pyo)
        end
    end
end

"""
    @pyattr class method

Macro for creating a method signature for a Julia Type
In Python world method(T::class) is analagous to class.method()
"""
macro pyattr(class, method)
    pyattr(class, method)
end

"""
    @pyattr class method orig_method

Macro for creating a method signature for a Julia Type
In Python world method(T::class) is analagous to class.method()
This will create a new method name which calls the orig_method
"""
macro pyattr(class, method, orig_method)
    pyattr(class, method, orig_method)
end

"""
    pyattr_set(types, methods...)

For each Julia type `T<:EEWrapped` in `types` and each method `m` in `methods`,
define a new function `m(t::T, args...)` that delegates to the underlying
pyobject wrapped by `t`.
"""
function pyattr_set(classes, methods...)
    for class in classes
        for method in methods
            @eval @pyattr($class, $method)
        end
    end
end

const ee_exports = OrderedSet()

"""
    Initialize(args...; kwargs...)

Function to initialize an EarthEngine session (analagous to ee.Initialize() from
the Python API). Accepts arguments and keywords from the Python ee.Initialize()
function. This function also dynamically builds the EE API and creates the methods
with signatures for each EE Type.
"""
function Initialize(args...; kwargs...)
    try
        copy!(ee, pyimport("ee"))
    catch err
        error(
            "The `earthengine-api` package could not be imported. You must install the Python earthengine-api before using this package. The error was $err",
        )
    end

    try
        ee.Initialize(args...; kwargs...)

    catch
        @warn "Could not initialize an `ee` session. Trying authentication workflow..."

        try
            Authenticate()
            ee.Initialize()
        catch
            error(
                "Could not initialize an `ee` session or run authentication workflow. Please authenticate manually using the earthengine-api CLI (i.e. `\$ earthengine authenticate`",
            )
        end
    end

    # pull in the types and dynamically wrap things after initialization
    include("$(module_dir)/eetypes.jl")
    include("$(module_dir)/eefuncs.jl")

    for f in Symbol.(ee_exports)
        @eval export $f
    end
end

"""
    Authenticate()

Function to execute the EarthEngine authetication workflow (analgous to
ee.Authenticate() in the Python API). This function should only be executed
once if the EE API has not be used before.
"""
function Authenticate(args...; kwargs...)
    try
        ee.Authenticate(args...; kwargs...)
    catch err
        @error "Could not authenticate `ee`...see following error: $err"
    end
end

# create an abreviated variable for the module
const EE = EarthEngine

export EE, ee, Initialize

end # module


text:
#    Copyright 2015 Raytheon BBN Technologies
#
#     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.

module SchattenNorms

export snorm, nucnorm, trnorm, specnorm, fnorm, dnorm

"""
Computes the nuclear norm of a matrix `m`.
"""
function nucnorm(m::AbstractMatrix)
    norm(svdvals(m),1)
end

"""
Computes the trace norm of a matrix `m`.
"""
trnorm(m::AbstractMatrix) = nucnorm(m)

"""
Computes the Frobenius norm of a matrix `m`.
"""
fnorm(m::AbstractMatrix) = vecnorm(m,2)

"""
Computes the spectral norm of a matrix `m` (i.e., the maximum singular value).
"""
function specnorm(m::AbstractMatrix)
    return norm(svdvals(m),Inf)
end

"""
Computes the `p`-Schatten norm of a matrix `m`.
"""
function snorm(m::AbstractMatrix,p=2.0)
    p == 2.0 && return vecnorm(m,2)
    return norm(svdvals(m),p)
end

include("dnorm.jl")

end # module


text:
function weibull(x, a, b)
    y = (a/b) .* (x./b).^(a-1) .* exp.(.-(x./b).^a)
    ys  = y ./ maximum(y)
    return y, ys
end

function carthesian2spherical(x::T, y::T, z::T) where {T<:AbstractFloat}
    theta = atan(sqrt(x*x + y*y), z)
    phi = atan(x, y)
    if phi < 0.0
        phi += 2.0*pi
    end
    return theta, phi
end

function spherical2carthesian(theta::T, phi::T) where {T<:AbstractFloat}
    x = sin(theta) * sin(phi)
    y = sin(theta) * cos(phi)
    z = cos(theta)
    return x, y, z
end

function relazi(vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, sx::T, sy::T, sz::T) where {T<:AbstractFloat}
    # v=viewVector; n=normalVector; h=halfVector...
    ux, uy, uz = cross_scalars(nx, ny, nz, sx, sy, sz)
    vx, vy, vz = cross_scalars(nx, ny, nz, vx, vy, vz)
    nsx, nsy, nsz = normalize_scalars(ux, uy, uz)
    nvx, nvy, nvz = normalize_scalars(vx, vy, vz)
    dotprod = nsx*nvx + nsy*nvy + nsz*nvz
    if dotprod > 1.0 # it would only be slightly higher by machine precision, e.g. 1.0000000000000002
        dotprod = 1.0
    end
    if dotprod < 0.0 # similar to the case above...
        dotprod = 0.0
    end
    a = acos(dotprod)
    return a
end

function translate(ray::Ray{T}, xoff::T, yoff::T, zoff::T) where {T<:AbstractFloat}
    x = ray.x + xoff
    y = ray.y + yoff
    z = ray.z + zoff
    tray = Ray(x, y, z, ray.dx, ray.dy, ray.dz)
    return tray
end

function transform(ray::Ray{T}, xoff::T, yoff::T, zoff::T, beta::T, gamma::T) where {T<:AbstractFloat}
    # to transform from world space into object space
    # (values for xoff, yoff, ..., gamma, should be of opposite sign to how they are defined in the scene specification file)

    # the translation
    x = ray.x + xoff
    y = ray.y + yoff
    z = ray.z + zoff

    # rotate about z-axis
    x_ = x * cos(gamma) + y * -sin(gamma)
    y  = x * sin(gamma) + y * cos(gamma)
    dx_ = ray.dx * cos(gamma) + ray.dy * -sin(gamma)
    dy  = ray.dx * sin(gamma) + ray.dy * cos(gamma)

    # rotate about y-axis
    x = x_ * cos(beta) + z * sin(beta)
    z = x_ * -sin(beta) + z * cos(beta)
    dx = dx_ * cos(beta) + ray.dz * sin(beta)
    dz = dx_ * -sin(beta) + ray.dz * cos(beta)

    tray = Ray(x, y, z, dx, dy, dz)
    return tray
end

function transform(x::T, y::T, z::T, xoff::T, yoff::T, zoff::T, beta::T, gamma::T) where {T<:AbstractFloat}
    # to transform x,y,z from object space into worldspace
    # rotate about y-axis
    t = x * cos(beta) + z * sin(beta)
    c = x * -sin(beta) + z * cos(beta)

    # rotate about z-axis
    a = t * cos(gamma) + y * -sin(gamma)
    b = t * sin(gamma) + y * cos(gamma)

    # and the translation
    a = a + xoff
    b = b + yoff
    c = c + zoff
    return a, b, c
end

# functions using Hook

function translate_forward(x::T, y::T, z::T, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    x = x + hook.xoff
    y = y + hook.yoff
    z = z + hook.zoff
    return x, y, z
end

function translate_inverse(x::T, y::T, z::T, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    x = x - hook.xoff
    y = y - hook.yoff
    z = z - hook.zoff
    return x, y, z
end

function translate_forward(x::AbstractArray{T}, y::AbstractArray{T}, z::AbstractArray{T}, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    x = x .+ hook.xoff
    y = y .+ hook.yoff
    z = z .+ hook.zoff
    return x, y, z
end

function translate_inverse(x::AbstractArray{T}, y::AbstractArray{T}, z::AbstractArray{T}, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    x = x .- hook.xoff
    y = y .- hook.yoff
    z = z .- hook.zoff
    return x, y, z
end

function rotate_forward(x::T, y::T, z::T, hook::Hook{T,S}) where {T <: AbstractFloat, S<:Integer}
    # rotate about y-axis
    nx_ = x * cos(hook.beta) + z * sin(hook.beta)
    nz =  x * -1.0 * sin(hook.beta) + z * cos(hook.beta)

    # rotate about z-axis
    nx =  nx_ * cos(hook.gamma) + y * -1.0 * sin(hook.gamma)
    ny =  nx_ * sin(hook.gamma) + y * cos(hook.gamma)
    return nx, ny, nz
end

function rotate_inverse(x::T, y::T, z::T, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    inv_beta = -1.0 * hook.beta
    inv_gamma = -1.0 * hook.gamma

    # rotate about z-axis (to get segment tip on y=0 plane)
    nx_1 = x * cos(inv_gamma) + y * -1.0 * sin(inv_gamma)
    ny = x * sin(inv_gamma) + y * cos(inv_gamma)

    # rotate about y-axis (to get segment tip on z=0 plane)
    nx = nx_1 * cos(inv_beta) + z * sin(inv_beta)
    nz = nx_1 * -1.0 * sin(inv_beta) + z * cos(inv_beta)
    return nx, ny, nz
end

function rotate_forward(x::AbstractArray{T}, y::AbstractArray{T}, z::AbstractArray{T}, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    # rotate about y-axis
    nx_ = x .* cos(hook.beta) .+ z .* sin(hook.beta)
    nz =  x .* -1.0 .* sin(hook.beta) .+ z .* cos(hook.beta)

    # rotate about z-axis
    nx =  nx_ .* cos(hook.gamma) .+ y .* -1.0 .* sin(hook.gamma)
    ny =  nx_ .* sin(hook.gamma) .+ y .* cos(hook.gamma)
    return nx, ny, nz
end

function rotate_inverse(x::AbstractArray{T}, y::AbstractArray{T}, z::AbstractArray{T}, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    inv_beta = -1.0 * hook.beta
    inv_gamma = -1.0 * hook.gamma

    # rotate about z-axis (to get segment tip on y=0 plane)
    nx_1 = x .* cos(inv_gamma) .+ y .* -1.0 .* sin(inv_gamma)
    ny = x .* sin(inv_gamma) .+ y .* cos(inv_gamma)

    # rotate about y-axis (to get segment tip on z=0 plane)
    nx = nx_1 .* cos(inv_beta) .+ z .* sin(inv_beta)
    nz = nx_1 .* -1.0 .* sin(inv_beta) .+ z .* cos(inv_beta)
    return nx, ny, nz
end

function transform_forward(x::T, y::T, z::T, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    x, y, z = rotate_forward(x, y, z, hook)
    x, y, z = translate_forward(x, y, z, hook)
    return x, y, z
end

function transform_inverse(x::T, y::T, z::T, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    x, y, z = translate_inverse(x, y, z, hook)
    x, y, z = rotate_inverse(x, y, z, hook)
    return x, y, z
end

function transform_forward(x::AbstractArray{T}, y::AbstractArray{T}, z::AbstractArray{T}, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    x, y, z = rotate_forward(x, y, z, hook)
    x, y, z = translate_forward(x, y, z, hook)
    return x, y, z
end

function transform_inverse(x::AbstractArray{T}, y::AbstractArray{T}, z::AbstractArray{T}, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    x, y, z = translate_inverse(x, y, z, hook)
    x, y, z = rotate_inverse(x, y, z, hook)
    return x, y, z
end

function transform_forward(ray::Ray{T}, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    x, y, z = rotate_forward(ray.x, ray.y, ray.z, hook)
    x, y, z = translate_forward(x, y, z, hook)
    dx, dy, dz = rotate_forward(ray.dx, ray.dy, ray.dz, hook)
    transRay::Ray = Ray(x, y, z, dx, dy, dz)
    return transRay
end

function transform_inverse(ray::Ray{T}, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    x, y, z = translate_inverse(ray.x, ray.y, ray.z, hook)
    x, y, z = rotate_inverse(x, y, z, hook)
    dx, dy, dz = rotate_inverse(ray.dx, ray.dy, ray.dz, hook)
    transRay::Ray = Ray(x, y, z, dx, dy, dz)
    return transRay
end

# functions using transformation matrix

function transform(x::T, y::T, z::T, matrix::AbstractArray{T,2}) where {T<:AbstractFloat}
    xnew = x*matrix[1,1] + y*matrix[1,2] + z*matrix[1,3] + matrix[1,4]
    ynew = x*matrix[2,1] + y*matrix[2,2] + z*matrix[2,3] + matrix[2,4]
    znew = x*matrix[3,1] + y*matrix[3,2] + z*matrix[3,3] + matrix[3,4]
    return xnew, ynew, znew
end

function transform(ray::Ray{T}, matrix::AbstractArray{T,2}) where {T<:AbstractFloat, S<:Integer}
    x2 = ray.x + ray.dx
    y2 = ray.y + ray.dy
    z2 = ray.z + ray.dz
    x, y, z = transform(ray.x, ray.y, ray.z, matrix)
    u, v, w = transform(x2, y2, z2, matrix)
    dx = u - x
    dy = v - y
    dz = w - z
    #dx, dy, dz = normalize(dx, dy, dz)   # DON'T APPLY NORMALIZE !!!
    transformedRay::Ray = Ray(x, y, z, dx, dy, dz)
    return transformedRay
end

# function using Transform (<:Asset) object

function transform_forward(x::T, y::T, z::T, t::Transform{T,S}) where {T<:AbstractFloat, S<:Integer}
    return transform(x, y, z, t.matrix)
end

function transform_inverse(x::T, y::T, z::T, t::Transform{T,S}) where {T<:AbstractFloat, S<:Integer}
    return transform(x, y, z, t.inverse)
end

function transform_forward(ray::Ray{T}, t::Transform{T,S}) where {T<:AbstractFloat, S<:Integer}
    return transform(ray, t.matrix)
end

function transform_inverse(ray::Ray{T}, t::Transform{T,S}) where {T<:AbstractFloat, S<:Integer}
    return transform(ray, t.inverse)
end

max_stars_repo_name: martinvanleeuwen/RenderJay

max_issues_repo_name: martinvanleeuwen/RenderJay

max_forks_repo_name: martinvanleeuwen/RenderJay


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to compute the printable width of the string.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export printable_textwidth, printable_textwidth_per_line

"""
    printable_textwidth(str::AbstractString)

Return the text width of `str` considering only the printable characters, *i.e.*
removing all ANSI escape sequences related to decorations.

!!! note
    Characters like `\\n` and `\\t` are treated as normal characters.
"""
function printable_textwidth(str::AbstractString)
    return str |> remove_decorations |> textwidth
end

"""
    printable_textwidth_per_line(str::AbstractString)

Return a vector with the printable textwidth of each line in `str`. The lines
are split considering the character `\n`.
"""
function printable_textwidth_per_line(str::AbstractString)
    lines = split(str, '\n')
    num_lines = length(lines)
    lines_width = zeros(Int, num_lines)

    @inbounds for k in 1:num_lines
        lines_width[k] = printable_textwidth(lines[k])
    end

    return lines_width
end


text:
module Statifiers

using Random
using StatsBase: std, skewness, kurtosis, variation, sem, mad
using StatsBase: entropy, summarystats, autocor, pacf, rle, quantile
using Dates
using DataFrames
using Statistics

using ..AbsTypes
using ..Utils
import ..AbsTypes: fit, fit!, transform, transform!

export fit, fit!,transform, transform!
export Statifier,tsmlfullstat

# Transforms instances with nominal features into one-hot form
# and coerces the instance matrix to be of element type Float64.

"""
    Statifier(Dict(
       :processmissing => true
    ))

Outputs summary statistics such as mean, median, quartile, entropy, kurtosis, skewness, etc.
with parameter:

- `:processmissing` => `boolean` to indicate whether to include `missing` data stats.

Example:

    dt=[missing;rand(1:10,3);missing;missing;missing;rand(1:5,3)]
    dat = DataFrame(Date= DateTime(2017,12,31,1):Dates.Hour(1):DateTime(2017,12,31,10) |> collect,
                    Value = dt)

    statfier = Statifier(Dict(:processmissing=>false))

    fit!(statfier,dat)
    results=transform!(statfier,dat)

Implements: `fit!`, `transform!`
"""
mutable struct Statifier <: Transformer
   name::String
   model::Dict{Symbol,Any}

   function Statifier(args=Dict())
      default_args = Dict(
          :name => "stfr",
          :processmissing => true
      )
      cargs=nested_dict_merge(default_args,args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      new(cargs[:name],cargs)
   end
end

"""
    fit!(st::Statifier, features::T=[], labels::Vector=[])

Validate argument to make sure it's a 2-column format.
"""
function fit!(st::Statifier, features::DataFrame=DataFrame(), labels::Vector=[])::Nothing
   typeof(features) <: DataFrame || throw(ArgumentError("Statifier.fit!: data should be a dataframe: Date,Val "))
   ncol(features) == 2 || throw(ArgumentError("dataframe must have 2 columns: Date, Val"))
   return nothing
end

function fit(st::Statifier, features::DataFrame=DataFrame(), labels::Vector=[])::Statifier
   fit!(st,features,labels)
   return deepcopy(st)
end

"""
    transform!(st::Statifier, features::T=[])

Compute statistics.
"""
function transform!(st::Statifier, features::DataFrame=DataFrame())::DataFrame
   features != [] || return DataFrame()
   typeof(features) <: DataFrame || throw(ArgumentError("Statifier.fit!: data should be a dataframe: Date,Val "))
   ncol(features) == 2 || throw(ArgumentError("dataframe must have 2 columns: Date, Val"))
   sum(names(features) .== ("Date","Value"))  == 2 || throw(ArgumentError("wrong column names"))
   fstat = tsmlfullstat(features.Value)
   timestat = timevalstat(features)
   if st.model[:processmissing] == true
      # full namedtuple: stat1,stat2,bstat
      hcat(timestat,fstat...)
   else
      hcat(timestat,fstat.stat1,fstat.stat2)
   end
end

function transform(st::Statifier, features::DataFrame=DataFrame())::DataFrame
   return transform!(st,features)
end

function timevalstat(features::DataFrame)
   ldates=features.Date |> skipmissing |> collect
   totalhours = sum(diff(ldates)).value/1000/3600 # to hours
   lcount = length(ldates)
   timestart=first(features.Date)
   timeend=last(features.Date)
   sfreq = totalhours/lcount
   dftime = DataFrame(tstart=timestart,tend=timeend,sfreq=sfreq)
end

function tsmlfullstat(dat::Vector)
   data = skipmissing(dat) |> collect
   lcount = length(data); lmax = maximum(data); lmin = minimum(data)
   lsm = summarystats(data)
   q1 = quantile(data,0.1); q2 = quantile(data,0.2)
   q8 = quantile(data,0.8); q9 = quantile(data,0.9)
   lks = kurtosis(data); lsk = skewness(data); lvar = variation(data)
   lsem = sem(data)
   lentropy = entropy(data)
   # assume 24-hour period
   _autolags = 1:minimum([24,length(data)-1])
   _pacflags = 1:minimum([24,div(length(data),2)-1])
   lautocor = autocor(data,_autolags) .|> abs2 |> sum |> sqrt
   lpacf = pacf(data,_pacflags;method=:yulewalker) .|> abs2 |> sum |> sqrt
   lbrle = rlestatmissingblocks(dat)
   df1=DataFrame(count=lcount,max=lmax,min=lmin,median=lsm.median,
                 mean=lsm.mean,q1=q1,q2=q2,q25=lsm.q25,
                 q75=lsm.q75,q8=q8,q9=q9)
   df2=DataFrame(kurtosis=lks,skewness=lsk,variation=lvar,entropy=lentropy,
                 autocor=lautocor,pacf=lpacf)
   df3=DataFrame(bmedian=lbrle.bmedian,bmean=lbrle.bmean,bq25=lbrle.bq25,bq75=lbrle.bq75,
                 bmin=lbrle.bmin,bmax=lbrle.bmax)
   return (stat1=df1,stat2=df2,bstat=df3)
end

function rlestatmissingblocks(dat::Vector{T}) where {T<:Union{AbstractFloat,Integer,Missing}}
   data = deepcopy(dat)
   # replace missing with dummy value for rle processing
   indxmissing=findall(x->ismissing(x),data)
   dummy = -99999
   if eltype(data) <: Union{AbstractFloat,Missing}
      data[indxmissing] .= dummy
   elseif eltype(data) <: Union{Integer,Missing}
      data[indxmissing] .= dummy
   else
      error("not a valid type")
   end
   _rle = rle(data)
   indx=findall(x->x == dummy,_rle[1])
   blocks = _rle[2][indx]
   _sm = summarystats(blocks)
   return (bq25=_sm.q25,bmean=_sm.mean,bmedian=_sm.median,bq75=_sm.q75,bmax=_sm.max,bmin=_sm.min,bmiss=_sm.nmiss)
end

end