text:
# Tue 23 Oct 2012 07:10:59 PM EDT
#
# NeHe Tut 4 - Rotate colored (rainbow) triangle and colored (cyan) square
#
# Q - quit

# load necessary GL/SDL routines

using OpenGL
@OpenGL.version "1.0"
@OpenGL.load
using SDL

# initialize variables

bpp                = 16
wintitle           = "NeHe Tut 4"
icontitle          = "NeHe Tut 4"
width              = 640
height             = 480

rtri               = 0.0
rquad              = 0.0

keystate_checked   = false
lastkeycheckTime   = 0
key_repeatinterval = 75 #ms

# open SDL window with an OpenGL context

SDL_Init(SDL_INIT_VIDEO)
#videoInfo = SDL_GetVideoInfo()
videoFlags = (SDL_OPENGL | SDL_GL_DOUBLEBUFFER | SDL_HWPALETTE | SDL_RESIZABLE)
#if videoInfo.hw_available
    videoFlags = (videoFlags | SDL_HWSURFACE)
#else
    #videoFlags = (videoFlags | SDL_SWSURFACE)
#end
#if videoInfo.blit_hw
    videoFlags = (videoFlags | SDL_HWACCEL)
#end
SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER, 1)
SDL_SetVideoMode(width, height, bpp, videoFlags)
SDL_WM_SetCaption(wintitle, icontitle)

glViewport(0, 0, width, height)
glClearColor(0.0, 0.0, 0.0, 0.0)
glClearDepth(1.0)
glDepthFunc(GL_LESS)
glEnable(GL_DEPTH_TEST)
glShadeModel(GL_SMOOTH)

glMatrixMode(GL_PROJECTION)
glLoadIdentity()

gluPerspective(45.0,width/height,0.1,100.0)

glMatrixMode(GL_MODELVIEW)

# main drawing loop

while true
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glLoadIdentity()

    glTranslate(-1.5,0.0,-6.0)
    glRotate(rtri,0.0,1.0,0.0)

    glBegin(GL_POLYGON)
      glColor(1.0,0,0)
      glVertex(0.0,1.0,0.0)
      glColor(0,1.0,0)
      glVertex(1.0,-1.0,0.0)
      glColor(0,0,1.0)
      glVertex(-1.0,-1.0,0.0)
    glEnd()

    glLoadIdentity()

    glTranslate(1.5,0.0,-6.0)
    glRotate(rquad,1.0,0.0,0.0)

    glTranslate(0.8,0,0)

    glColor(0.5,0.5,1.0)
    glBegin(GL_QUADS)
        glVertex(-1.0,1.0,0.0)
        glVertex(1.0,1.0,0.0)
        glVertex(1.0,-1.0,0.0)
        glVertex(-1.0,-1.0,0.0)
    glEnd()

    rtri  +=0.2
    rquad -=0.2

    SDL_GL_SwapBuffers()

    SDL_PumpEvents()
    if SDL_GetTicks() - lastkeycheckTime >= key_repeatinterval
        keystate         = SDL_GetKeyState()
        keystate_checked = true
        lastkeycheckTime = SDL_GetTicks()
    end

    # Sampling rates for event processing are so fast that a single key press
    # lasts through several iterations of this loop.  This means that one press
    # can be seen as 50 or more presses by the SDL event system.  To correct
    # for this, we only check keypresses every 75ms.

    if keystate_checked == true
        if keystate[SDLK_q] == true
            break
        end
        keystate_checked = false
    end
end


text:
module SequenceMutationDistance

    include("include_all.jl")

    export

    #smdfuncs
    dist_matrix,
    dist_matrix_mt,
    symmetric_dist_matrix,
    symmetric_dist_matrix_mt,
    smd_mf,
    smd_sum,
    smd_distance_wrapper,
    fst_wrapper,
    permutation_test

end  # module

ext: jl


text:
"""
    fba(model,[direction = "max", objective = 0)

Solves the LP : max/min c'*v
                  subject to S*v = b
                  lb <= v <= ub

### Examples
    fbasolution = fba(model)

to maximize the default objective

    fbasolution = fba(model, objective="ACALD")

to maximize the "ACALD" reaction

    fbasolution = fba(model, objective=4)

to maximize reaction 4

Return a FBAsolution object which has the following fields:
* obj - Objective value
* v - Primal solution
* slack - Dual solution
* rcost - Reduced costs
* success - true if solution was a success, false otherwise
* info - Solver info
"""
function fba(model::Model; direction::String = "max", objective = 0)
    if typeof(objective) <: String
        if any(model.rxns .== objective)
            objective = first(find(model.rxns .== objective))
        else
            warn("Objective reaction not found in model.rxns")
            return
        end
    end

    fba(model, direction, objective)
end

function fba(model::Model, direction::String = "max", objective = 0)
	lp = setup_lp(model, direction);

    if objective != 0
        if objective > length(model.c)
            warn("Objective is larger than number of reactions")
            return
        end
        change_objective_coef(lp, objective)
    end
	solve(lp)

	return construct_lp_solution(lp)
end

"""
    fva(model,[optPercentage = 100])

Determine the flux variability of the model

### Optional arguements:
    optPercentage: 0 - 1 (fraction)

Used to fix the biomass production rate anywhere between 0 - 1

#### Examples
Run flux variability with biomass lower bound at 90% of maximum

    min, max, minX, maxX = fva(model, 0.9)

Return a min flux vector, max flux vector, minimum primal flux vector and maximum
primal flux vector
"""
function fva(model::Model; optPercentage::Number = 1, flux_matrix::Bool = false)
	fva(model, optPercentage, flux_matrix)
end

function fva(model::Model, optPercentage::Number = 1, flux_matrix::Bool = false)

    if length(procs()) != 1
		return PCBM.fast_fva(model, optPercentage, flux_matrix)
	end

	num_rxns = model.S.n
	objective_index = find(model.c)
	progress_meter = ProgressMeter.Progress(num_rxns*2, 0.1, "\n\u1b[1Ffva:", 50)

	lp = setup_lp(model, "max")

	# fix the biomass lower bound according to "optPercentage"
	solve(lp)
	new_lb = get_solution(lp)[objective_index]
	new_lb *= optPercentage

	map((x,y) -> set_col_bounds(lp, x, y, model.ub[x]), objective_index, new_lb)
	map(x -> set_objective(lp, x, 0.0), objective_index)

	set_direction(lp, "min")
	minFlux, minV =  variability_runner(lp, num_rxns, progress_meter, flux_matrix)

	set_direction(lp, "max")
	maxFlux, maxV =  variability_runner(lp, num_rxns, progress_meter, flux_matrix)

    if flux_matrix
       return minFlux,maxFlux, minV, maxV
    else
	   return minFlux,maxFlux
    end
end

function variability_runner(lp, num_rxns, progress_meter, flux_matrix)
	flux = Array(Float64, num_rxns)
	v = Array(Float64, num_rxns, num_rxns)
	for i = 1:num_rxns

		set_objective(lp,i)
		solve(lp)

		flux[i] = get_objective_value(lp)

        if flux_matrix
		  v[:,i] = get_solution(lp)
        end

		set_objective(lp,i, 0.0)
		next!(progress_meter)
	end

	return flux, v
end

# -------------------------------------------------------------------
# -------------------------------------------------------------------

function jump_fba(model)
  # Flux balance analysis
  #    sol=fba(model) performs flux balance Analysis
  #
  # Ref.  Savinell J, Palsson . J. Theor. Biol. 1992, 154(4):421-454.

  # THINK: Add side constraints, e.g. minimize 1 and 2 norms
  # 0-norm corresponds to min_set_rxns

  nrxns, nmets, _ = statistics(model)
  m, v, constr = create_fba_model(model)
  status = solve(m)
  info=SolverInfo(cobra_lp_solver, status)
  if status == :Optimal
    return FBAsol(getobjectivevalue(m), getvalue(v), getdual(constr), getdual(v), true, info)
  end
  return FBAsol(0, [], [], [], false, info)
end

function getmaxobjective(model)
  sol = fba(model)
  if !sol.success
    error("Unable to solve to optimality")
  end
  sol.obj
end

function fba_norm2(model::Model)
  # Flux balance analysis - minimum 2-norm solution

  optobj=getmaxobjective(model)
  nrxns, nmets, _ = statistics(model)
  m = createqpmodel()
  @variable(m, model.lb[i] <= v[i=1:nrxns] <= model.ub[i])
  @objective(m, Min, sum{v[i]^2, i=1:nrxns})
  @constraint(m, model.S*v .== model.m.b)
  @constraint(m, dot(v, model.c) >= optobj)
  status = solve(m)
  info = SolverInfo(cobra_qp_solver, status)
  if status == :Optimal
    return FBAsol(optobj, getvalue(v), [], [], true, info)
  end
  return FBAsol(0, [], [], [], false, info)
end

function fba_norm1(model::Model)
  # Flux balance analysis - minimum 1-norm solution

  # THINK:A lot of this is very similar to fba_norm2
  global cobra_warning_level
  optobj=getmaxobjective(model)
  nrxns, nmets, _ = statistics(model)
  m = createlpmodel()
  @variable(m, model.r.lb[i] <= v[i=1:nrxns] <= model.r.ub[i])
  @variable(m, z[i=1:nrxns])
  @objective(m, Min, sum(z))
  @constraint(m, model.s*v .== model.m.b)
  @constraint(m, dot(v, model.r.c) >= optobj)
  @constraint(m, z .>= v)
  @constraint(m, z .>= -v)
  status = solve(m)
  info = SolverInfo(cobra_lp_solver, status)
  if status == :Optimal
    return FBAsol(optobj, getvalue(v), [], [], true, info)
  else
    if cobra_warning_level
      warn("Unable to solve to optimality, status: $status")
    end
    return FBAsol(0, [], [], [], false, info)
  end
end

#function fba_norm0(model::Model, rxns)
function fba_norm0(model::Model)
  # Flux balance analysis - minimum number of nonzero fluxes (min 0-norm)

  # THINK: TBD Implement subset of reactions

  optobj=getmaxobjective(model)
  nrxns, nmets, _ = statistics(model)
  m = createmilpmodel()
  @variable(m, v[i=1:nrxns])
  @variable(m, y[i=1:nrxns], Bin) # THINK: subset
  @objective(m, Min, sum(y))
  @constraint(m, model.S*v .== model.b)
  @constraint(m, dot(v, model.c) >= optobj)
  @constraint(m, tieub[i=1:nrxns], v[i] <= model.ub[i]*y[i])
  @constraint(m, tielb[i=1:nrxns], model.lb[i]*y[i] <= v[i])
  status = solve(m)
  info=SolverInfo(cobra_mip_solver, status)
  if status == :Optimal # THINK: Suboptimal solutions?
    return FBAsol(optobj, getvalue(v), [], [], true, info)
  else
    if cobra_warning_level
      warn("Unable to solve to optimality, status: $status")
    end
    return FBAsol(0, [], [], [], false, info)
  end
end

#function fba_norm0(model::Model)=fba_norm0(model, collect(model.r.name))
function jump_fva(model::Model, percent=100)
  # Flux variability  analysis
  #    minFlux, maxFlux = fva(model, percent) performs flux variability analysis

  # Ref. Mahadevan R, Schilling C. Metab. Eng. 2003, 5(4):264-276.

  # THINK: Side constraints?

  global cobra_warning_level, cobra_tolerance
  nrxns, nmets, _ = statistics(model)
  m, x, constr = create_fba_model(model)
  status = solve(m)
  if status != :Optimal
    if cobra_warning_level
      warn("Unable to solve to optimality, status: $status")
    end
    return [], [], false, status
  end
  if percent > 0
    # Lower bound the cellular objective
    obj_lb = floor(getobjectivevalue(m) / cobra_tolerance) * cobra_tolerance * percent / 100.0;
    @constraint(m, dot(x, model.c) >= obj_lb)
  end

  minFlux, maxFlux = zeros(Float64, nrxns), zeros(Float64, nrxns)
  # Solve all maximization problems first, then the minimization problems
  for i=1:nrxns
    c = zeros(nrxns)
    c[i] = 1
    @objective(m, Max, dot(x, c))
    status = solve(m)
    maxFlux[i] = getobjectivevalue(m)
  end
  for i=1:nrxns
    c = zeros(nrxns)
    c[i] = 1
    @objective(m, Min, dot(x, c))
    status = solve(m)
    minFlux[i] = getobjectivevalue(m)
  end
  return minFlux, maxFlux, true, status
end

function moma(model_wt::Model, model_del::Model)
  # Minimization of metabolic adjustment (MOMA)
  #    moma(model_wt, model_del)

  # Ref: Segre et al. PNAS (2002), 99(23), 15112-15117

  # This formulation is identical to the one in the COBRA 2.0 toolbox. It
  # differs from Segre et al. in such a way that the wild-type
  # and mutant flux distributions are determined simultaneously. This is done
  # to avoid dependence on solver/solver settings which arises because of the
  # existence of alternative optima in the FBA LP.

  optobj=getmaxobjective(model)

  # Solve the min norm problem
  nrxns, nmets, _ = statistics(model) # THINK: Silly function
  m = createqpmodel()
  @variable(m, model_wt.r.lb[i] <= v_wt[i=1:nrxns] <= model_wt.r.ub[i])
  @variable(m, model_del.r.lb[i] <= v_del[i=1:nrxns] <= model_del.r.ub[i])
  @objective(m, Min, sum{(v_wt[i] - v_del[i])^2, i=1:nrxns})
  @constraint(m, model_wt.s*v_wt .== model_wt.m.b)
  @constraint(m, model_del.s*v_del .== model_del.m.b)
  @constraint(m, dot(model_wt.r.c, v_wt) == optobj)
  status = solve(m)
  info=SolverInfo(cobra_qp_solver, status)
  return MOMAsol(getobjectivevalue(m), getvalue(v_wt), getvalue(v_del), (status == :Optimal), info)
end

function create_fba_model(model)
  # Create FBA linear program
  nrxns, nmets, _ = statistics(model)
  m = createlpmodel()
  @variable(m, model.lb[i] <= v[i=1:nrxns] <= model.ub[i])
  @objective(m, Max, dot(v,model.c))
  @constraint(m, constr, model.S*v .== model.b)
  solve(m)
  return m, v, constr
end

max_stars_repo_name: isebarn/CBM

max_issues_repo_name: isebarn/CBM

max_forks_repo_name: isebarn/CBM


text:
@info "Including exp.jl .."
include("../src/exp.jl")

function main()
    for d in [10,20,50,100],
        (prefix, model_fn,nsteps) in [
            ("EQ2", eq2_model_fn, 5e4),
            ("FC", ()->fc_model(d=d, ch=2, z=1024, nlayer=6), 3e4),
            ("CNN", ()->cnn_model(2), 3e4),
        ]

        specs = []
        for gtype in [:ER, :SF],
            k in [1]
            push!(specs, DataSpec(d=d, k=k, gtype=gtype,
                    noise=:Gaussian, mat=:CH3))
        end
        specs = Array{DataSpec}(specs)

        test_throttle = if prefix == "EQ2" 10 else 1 end

        @info "training .." prefix d
        expID = exp_train(specs, model_fn,
                          prefix="$prefix-ERSF-k1-d=$d", train_steps=nsteps,
                          test_throttle = test_throttle,
                          merge=true)
        @info "testing"
        exp_test(expID, specs, "TEST-gtype=$gtype")
    end
end

@info "Running main() .."
main()

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
using NLopt, ForwardDiff

x0 = zeros(10)
gcfg = ForwardDiff.GradientConfig(nothing,x0)
g!(grad, x) = ForwardDiff.gradient!(grad, integral_function, x, gcfg)
lbs  = -5ones(10)
lbs[6] = 0.0
ubs = 20ones(10)
ubs[6] = 1.0

function integral_function(x)
  a2, a3, b13, b22, b23, A021, A031, B021, B031, B032 = x
  c01 = 0
  c03 = 1
  c11 = 0
  c13 = 1
  A032 = 1 - A031
  b12 = (b22*(1 - c11 + b13*c11 - b13*c13))/(-1 + b23*c11 - b23*c13)
  b11 = 1 - b12 - b13
  b21 = - b22 - b23
  c12 = (-1 - b21*c11 - b23*c13)/b22
  a1 = 1 - a2 - a3
  c02 = A021
  G(A021,A031,A032,a1,a2,a3)
end

iter = 0
g_res = zeros(3)

function f(x::Vector, grad::Vector)
  if length(grad) > 0
    g!(grad,x)
  end
  integral = integral_function(x)
  global iter
  iter::Int += 1
  if mod(iter,1000) == 0
    global g_res
    g_res[1] = a2*B021 + a3*B031 + a3*B032 - 1
    g_res[2] = a2*B021^2 + a3*(B031 + B032)^2 - 3/2
    g_res[3] = A021*a2 + A031*a3 + A032*a3 - 1/2
    @show x
    @show g_res
    @show integral
  end
  integral
end

function G(A021,A031,A032,a1,a2,a3)
  tmps = zeros(Threads.nthreads())
  xmin = 6
  ymin = 4
  xs = -xmin:0.01:1
  ys = -ymin:0.01:4
  Threads.@threads for x in xs
    for y in ys
      z = x + im*y
      tmps[Threads.threadid()] += abs(1 + z*(a1 + a2 + A021*z*a2 + a3 + A032*z*a3 + (A031*z + A021*A032*z^2)*a3))<1
    end
  end
  tmp = sum(tmps)
  (tmp/(length(xs)*length(ys)))*((xmin+1)*(2ymin))
end

function constraints!(res,x,grad)
  a2, a3, b13, b22, b23, A021, A031, B021, B031, B032 = x
  c01 = 0
  c03 = 1
  c11 = 0
  c13 = 1
  A032 = 1 - A031
  b12 = (b22*(1 - c11 + b13*c11 - b13*c13))/(-1 + b23*c11 - b23*c13)
  b11 = 1 - b12 - b13
  b21 = - b22 - b23
  c12 = (-1 - b21*c11 - b23*c13)/b22
  a1 = 1 - a2 - a3
  c02 = A021
  if length(grad)>0
    grad .= 0.0
    grad[1,1] = B021
    grad[2,1] = B031 + B032
    grad[12,1] = a2
    grad[13,1] = a3
    grad[14,1] = a3
    grad[1,2] = B021^2
    grad[2,2] = (B031 + B032)^2
    grad[12,2] = 2a2*B021
    grad[13,2] = 2a3*(B031+B032)
    grad[14,2] = 2a3*(B031+B032)
    grad[1,3] = A021
    grad[2,3] = A031 + A032
    grad[9,3] = a2
    grad[10,3] = a3
    grad[11,3] = a3
  end
  res[1] = a2*B021 + a3*B031 + a3*B032 - 1
  res[2] = a2*B021^2 + a3*(B031 + B032)^2 - 3/2
  res[3] = A021*a2 + A031*a3 + A032*a3 - 1/2
end

function c_constraints!(res,x,grad)
  a2, a3, b13, b22, b23, A021, A031, B021, B031, B032 = x
  c01 = 0
  c03 = 1
  c11 = 0
  c13 = 1
  A032 = 1 - A031
  b12 = (b22*(1 - c11 + b13*c11 - b13*c13))/(-1 + b23*c11 - b23*c13)
  b11 = 1 - b12 - b13
  b21 = - b22 - b23
  c12 = (-1 - b21*c11 - b23*c13)/b22
  a1 = 1 - a2 - a3
  c02 = A021
  res[1] = c12 - 1
  res[2] = -c12
end

a2 = 1/6; a3 = 2/3; c01 = 0; c11 = 1; c13 = 0; b13 = 0;
b22 = 1; b23 = 0; A021 = 1; A031 = 1/4; A032 = 1/4;
B021 = 0; B031 = 1; B032 = 1/2
#x0 = [a2, a3, c01, c11, c13, b13, b22, b23, A021, A031, A032, B021, B031, B032] + rand(14)
x0 = rand(10)
x0[6] = max.(0,x0[6])
x0[6] = min.(1,x0[6])
opt = Opt(:LN_COBYLA, 10) #:LD_SLSQP, :LN_COBYLA (semi), :GN_ISRES, :LN_AUGLAG_EQ support equality constraints
lower_bounds!(opt, lbs)
upper_bounds!(opt, ubs)
max_objective!(opt, f)
equality_constraint!(opt::Opt,constraints!, [1e-12 for i in 1:3])
inequality_constraint!(opt::Opt,c_constraints!, [1e-12 for i in 1:2])
ftol_abs!(opt,1e-12)
xtol_rel!(opt,1e-12)
maxeval!(opt::Opt, Int(1e6))
maxf,maxx,ret = optimize(opt,x0)
println(maxf)

#maxx_str = "x = [$(maxx[1]),$(maxx[2]),$(maxx[3]),$(maxx[4]),$(maxx[5]),$(maxx[6]),$(maxx[7]),$(maxx[8]),$(maxx[9]),$(maxx[10]),$(maxx[11]),$(maxx[12]),$(maxx[13]),$(maxx[14])]"
#println(maxx_str)

### Result

x = maxx
a2, a3, b13, b22, b23, A021, A031, B021, B031, B032 = x
c01 = 0
c03 = 1
c11 = 0
c13 = 1
A032 = 1 - A031
b12 = (b22*(1 - c11 + b13*c11 - b13*c13))/(-1 + b23*c11 - b23*c13)
b11 = 1 - b12 - b13
b21 = - b22 - b23
c12 = (-1 - b21*c11 - b23*c13)/b22
a1 = 1 - a2 - a3
c02 = A021
G(A021,A031,A032,a1,a2,a3)
g_res = zeros(3); g_res_grad = zeros(14,3); c_res=zeros(4)
constraints!(g_res, x, g_res_grad)
g_res
g_res_grad
c_constraints!(c_res, x, g_res_grad)
c_res

value_str = """
  a1 = $a1;
  a2 = $a2;
  a3 = $a3;
  c01 = $c01;
  c02 = $c02;
  c03 = $c03;
  c11 = $c11;
  c12 = $c12;
  c13 = $c13;
  b11 = $b11;
  b12 = $b12;
  b13 = $b13;
  b21 = $b21;
  b22 = $b22;
  b23 = $b23;
  A021 = $A021;
  A031 = $A031;
  A032 = $A032;
  B021 = $B021;
  B031 = $B031;
  B032 = $B032;
  """
println(value_str)
#=
##################
# Test Integral
G(1,13/56,-1/56,13/4,5/4,-7/2) # 12.063
G(1,1/4,1/4,1/6,1/6,2/3) # 9.0872

###################
# Test constraints

a2 = 1/6; a3 = 2/3; c01 = 0; c11 = 1; c13 = 0; b13 = 0;
b22 = 1; b23 = 0; A021 = 1; A031 = 1/4; A032 = 1/4;
B021 = 0; B031 = 1; B032 = 1/2
x = [a2, a3, c01, c11, c13, b13, b22, b23, A021, A031, A032, B021, B031, B032]

a2, a3, c01, c11, c13, b13, b22, b23, A021, A031, A032, B021, B031, B032 = x
b12 = (b22*(1 - c11 + b13*c11 - b13*c13))/(-1 + b23*c11 - b23*c13)
b11 = 1 - b12 - b13
b21 = - b22 - b23
c12 = (-1 - b21*c11 - b23*c13)/b22
a1 = 1 - a2 - a3
c02 = A021
c03 = A031 + A032

a1 ≈ 1/6
a2 ≈ 1/6
a3 == 2/3
c01 == 0
c02 == 1
c03 == 1/2
c11 == 1
c12 == 0
c13 == 0
A021 == 1
A031 == 1/4
A032 == 1/4
B021 == 0
B031 == 1
B032 == 1/2
b11 == 1
b12 == 0
b13 == 0
b21 == -1
b22 == 1
b23 == 0

g_res = zeros(3); g_res_grad = zeros(14,3); c_res=zeros(4)
constraints!(g_res, x, g_res_grad)
c_constraints!(c_res, x, g_res_grad)
g_res
g_res_grad
c_res
=#


text:
using UnicodePlots # for better plotting, see Plots.jl
using GarbageMath # a joke package based on https://xkcd.com/2295/

# creates a set of "N" linearly spaced points on [-1,1]
N = 1000
x = LinRange(-1,1,N)

# creates a function
f(x) = exp(2*sin(pi*x))

# "broadcasts" f(x) so that it evaluates f(x[i]) for each x[i] ∈ x
fx = f.(x) # try typing "f(x)" into the REPL. Does it work?
display(lineplot(x,fx))

# testing the joke package
@show Precise(1.2) + Precise(3.5) # "@show" displays the output in the REPL
@show Precise(1) + Garbage(2)
@show Garbage(6)^2
@show sqrt(Garbage(5))


text:
using GeometricIntegrators.Solvers
using Test

A = [[+4.  +5.  -2.]
     [+7.  -1.  +2.]
     [+3.  +1.  +4.]]
b = [-14., +42., +28.]
x = [+4., -4., +5.]

function test_lu_solver(solver, A, b, x)
    # for T in (Float64, ComplexF64, Float32, ComplexF32) # TODO
    for T in (Float64, ComplexF64)
        AT = convert(Array{T,2}, A)
        bT = convert(Array{T,1}, b)
        xT = convert(Array{T,1}, x)

        lu = solver(AT, bT)
        factorize!(lu)
        solve!(lu)
        @test lu.b ≈ xT atol=1E-14
    end
end

test_lu_solver(LUSolverLAPACK, A, b, x)
test_lu_solver(LUSolver, A, b, x)


text:
## Global initialisation
using CompScienceMeshes
using BoundaryElements
import JLD

## Auxiliary scripts
include(Pkg.dir("CompScienceMeshes","examples","mesh_box_width_lid.jl"))
include(Pkg.dir("CompScienceMeshes","examples","matlab_patches.jl"))
include(Pkg.dir("BoundaryElements","examples","matlab_plotting.jl"))
@matlab cd("C:\\Temp")

## operators and excitations
const λ = 3.0
const κ = 2π/λ
const γ = im*κ
const μ = 1.0
const ϵ = 1.0
const η = sqrt(μ / ϵ)
const ω = κ / sqrt(μ * ϵ)
const α = 1 / γ

const I = BoundaryElements.Identity()
const N = BoundaryElements.NCross()
const T = BoundaryElements.MWSingleLayer3D(γ)
const K = BoundaryElements.MWDoubleLayer3D(γ)
const P = BoundaryElements.SingleLayerTrace(γ)
const nxK = BoundaryElements.DoubleLayerRotatedMW3D(γ)

const F = BoundaryElements.MWFarField3D(γ)
const Θ = linspace(0,π,100)
const Φ = [0.0]
const ffpts = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in Θ for ϕ in Φ]

const d = point(1,0,1); d = normalize(d)
const p = point(0,1,0); p = cross(p,d); p = normalize(p)
const E = BoundaryElements.planewavemw3d(direction = d, polarization = p, wavenumber = κ)
const H = 1 / (im*ω*μ) * curl(E)

const et = (n × E) × n
const hx = n × H
const ht = (n × H) × n

## define the solvers
function efie(m)
    X = raviartthomas(m)
    A = assemble(T, X, X)
    b = assemble(et, X)
    u = A \ b
    return u, X
end

function mxmfie(m)
    X = raviartthomas(m)
    Y = buffachristiansen(m)
    A = assemble(K + 0.5N, Y, X)
    b = assemble(ht, Y)
    u = A \ b
    return u, X
end

function clmfie(m)
    X = raviartthomas(m)
    A = assemble(nxK - 0.5I, X, X)
    b = assemble(hx, X)
    u = A \ b
    return u, X
end

##shorthand
farfield(u,X) = potential(F, ffpts, u, X)

## reference solution
m0 = meshcube(1.0, 0.04)
u0, X0 = efie(m0)
ff0 = farfield(u0, X0)
JLD.@save "refsolution.jld" u0

## load the reference solution
m0 = meshcube(1.0, 0.04)
X0 = raviartthomas(m0)
JLD.@load "refsolution.jld"
ff0 = farfield(u0, X0)

## the mfie solutions
hs = [0.2 0.15 0.14 0.13 0.12 0.11 0.1 0.08]
#hs = hs[1:2]

ff1, us1, Xs1 = [], [], []
ff2, us2, Xs2 = [], [], []
ff3, us3, Xs3 = [], [], []

for (i,h) in enumerate(hs)
    h = hs[i]

    m = meshcube(1.0, h)
    @show i, numcells(m)

    u1, X1 = efie(m)
    u2, X2 = clmfie(m)
    u3, X3 = mxmfie(m)

    push!(us1, u1); push!(Xs1, X1)
    push!(us2, u2); push!(Xs2, X2)
    push!(us3, u3); push!(Xs3, X3)

    push!(ff1, farfield(u1, X1))
    push!(ff2, farfield(u2, X2))
    push!(ff3, farfield(u3, X3))
end


text:
module GetPositionInstructions

  function get_position(direction, distance, verticle_distance, horizontal_distance)
    if direction == "forward"
      horizontal_distance += distance
    elseif direction == "up"
      verticle_distance -= distance
    else
      verticle_distance += distance
    end
    return verticle_distance, horizontal_distance
  end

  verticle_distance = 0
  horizontal_distance = 0

  instructions = readlines("input.txt", keep=true)
  function final_position(instructions, get_position, verticle_distance, horizontal_distance)
    for instruction in instructions
      direction, distance = split(instruction)
      distance = parse(Int64, distance)
      verticle_distance, horizontal_distance = get_position(direction, distance, verticle_distance, horizontal_distance)
    end
    return verticle_distance*horizontal_distance
  end

  print("The final position of the ship is: ", final_position(instructions, get_position, verticle_distance, horizontal_distance))
end
export final_position, get_position


text:
"""
    eulerivp(dudt,tspan,u0,n)

Apply Euler's method to solve the IVP u'=`dudt`(u,t) over the interval `tspan` with
u(`tspan[1]`)=`u0`, using `n` subintervals/steps. Return vectors of times and solution
values.
"""
function eulerivp(dudt,tspan,u0,n)
    a,b = tspan
    h = (b-a)/n
    t = [ a + i*h for i=0:n ]
    u = zeros(n+1)
    u[1] = u0
    for i = 1:n
      u[i+1] = u[i] + h*dudt(u[i],t[i])
    end
    return t,u
end

"""
    eulersys(dudt,tspan,u0,n)

Apply Euler's method to solve the vector-valued IVP u'=`dudt`(u,p,t) over the interval
`tspan` with u(`tspan[1]`)=`u0`, using `n` subintervals/steps.

(This is the same code as `eulerivp`. In the system case it returns a vector of vector
values as the solution.)
"""
function eulersys(dudt,tspan,u0,n)
    # Time discretization.
    a,b = tspan
    h = (b-a)/n
    t = [ a + i*h for i=0:n ]

    # Initial condition and output setup.
    u = Vector{typeof(u0)}(undef,n+1)
    u[1] = u0

    # The time stepping iteration.
    for i = 1:n
        u[i+1] = u[i] + h*dudt(u[i],t[i])
    end
    return t,u
end

"""
    ie2(dudt,tspan,u0,n)

Apply the Improved Euler method to solve the vector-valued IVP u'=`dudt`(u,p,t) over the
interval `tspan` with u(`tspan[1]`)=`u0`, using `n` subintervals/steps. Returns a vector
of times and a vector of solution values/vectors.
"""
function ie2(dudt,tspan,u0,n)
    # Time discretization.
    a,b = tspan
    h = (b-a)/n
    t = [ a + i*h for i=0:n ]

    # Initialize output.
    u = Vector{typeof(u0)}(undef,n+1)
    u[1] = u0

    # Time stepping.
    for i = 1:n
        uhalf = u[i] + h/2*dudt(u[i],t[i]);
        u[i+1] = u[i] + h*dudt(uhalf,t[i]+h/2);
    end
    return t,u
end

"""
    rk4(dudt,tspan,u0,n)

Apply "the" Runge-Kutta 4th order method to solve the vector-valued IVP u'=`dudt`(u,p,t)
over the interval `tspan` with u(`tspan[1]`)=`u0`, using `n` subintervals/steps.
Return a vector of times and a vector of solution values/vectors.
"""
function rk4(dudt,tspan,u0,n)
    # Time discretization.
    a,b = tspan
    h = (b-a)/n
    t = [ a + i*h for i=0:n ]

    # Initialize output.
    u = Vector{typeof(u0)}(undef,n+1)
    u[1] = u0

    # Time stepping.
    for i = 1:n
        k1 = h*dudt( u[i],      t[i]     )
        k2 = h*dudt( u[i]+k1/2, t[i]+h/2 )
        k3 = h*dudt( u[i]+k2/2, t[i]+h/2 )
        k4 = h*dudt( u[i]+k3,   t[i]+h   )
        u[i+1] = u[i] + (k1 + 2*(k2 + k3) + k4)/6
    end
    return t,u
end

"""
    rk23(dudt,tspan,u0,tol)

Apply adaptive embedded RK formula to solve the vector-valued IVP u'=`dudt`(u,p,t)
over the interval `tspan` with u(`tspan[1]`)=`u0`, with error tolerance `tol`.
Return a vector of times and a vector of solution values/vectors.
"""
function rk23(dudt,tspan,u0,tol)
    # Initialize for the first time step.
    t = tspan[1]
    u = [u0];   i = 1;
    h = 0.5*tol^(1/3)
    s1 = dudt(u0,t)

    # Time stepping.
    while t[i] < tspan[2]
        # Detect underflow of the step size.
        if t[i]+h == t[i]
            @warn "Stepsize too small near t=$(t[i])"
            break  # quit time stepping loop
        end

        # New RK stages.
        s2 = dudt( u[i]+(h/2)*s1,   t[i]+h/2   )
        s3 = dudt( u[i]+(3*h/4)*s2, t[i]+3*h/4 )
        unew2 = u[i] + h*(2*s1 + 3*s2 + 4*s3)/9   # 2rd order solution
        s4 = dudt( unew2, t[i]+h )
        err = h*(-5*s1/72 + s2/12 + s3/9 - s4/8)    # 2nd/3rd order difference
        E = norm(err,Inf)                           # error estimate
        maxerr = tol*(1 + norm(u[i],Inf))         # relative/absolute blend

        # Accept the proposed step?
        if E < maxerr     # yes
            t = [ t; t[i] + h]
            push!(u,unew2)
            i = i+1;
            s1 = s4;      # use FSAL property
        end

        # Adjust step size.
        q = 0.8*(maxerr/E)^(1/3)       # conservative optimal step factor
        q = min(q,4)                   # limit stepsize growth
        h = min(q*h,tspan[2]-t[i])     # don't step past the end
    end
    return t,u
end

"""
    ab4(dudt,tspan,u0,n)

Apply the Adams-Bashforth 4th order method to solve the vector-valued IVP u'=`dudt`(u,p,t)
over the interval `tspan` with u(`tspan[1]`)=`u0`, using `n` subintervals/steps.
"""
function ab4(dudt,tspan,u0,n)
    # Time discretization.
    a,b = tspan
    h = (b-a)/n
    t = [ a + i*h for i=0:n ]

    # Constants in the AB4 method.
    k = 4;    sigma = [55, -59, 37, -9]/24;

    # Find starting values by RK4.
    u = Vector{typeof(u0)}(undef,n+1)
    ts,us = rk4(dudt,[a,a+(k-1)*h],u0,k-1)
    u[1:k] = us[1:k]

    # Compute history of u' values, from newest to oldest.
    f = [ dudt(u[k-j],t[k-j]) for j=1:k-1  ]

    # Time stepping.
    for i = k:n
      f = [ dudt(u[i],t[i]), f[1:k-1]... ]   # new value of du/dt
      u[i+1] = u[i] + h*sum(f[j]*sigma[j] for j=1:k)       # advance one step
    end
    return t,u
end

"""
    am2(dudt,tspan,u0,n)

Apply the Adams-Moulton 2nd order method to solve the vector-valued IVP u'=`dudt`(u,p,t)
over the interval `tspan` with u(`tspan[1]`)=`u0`, using `n` subintervals/steps.
"""
function am2(dudt,tspan,u0,n)
    # Time discretization.
    a,b = tspan
    h = (b-a)/n
    t = [ a + i*h for i=0:n ]

    # Initialize output.
    u = Vector{typeof(u0)}(undef,n+1)
    u[1] = u0

    # Time stepping.
    for i = 1:n
        # Data that does not depend on the new value.
        known = u[i] + h/2*dudt(u[i],t[i])
        # Find a root for the new value.
        F = z -> z - h/2*dudt(z,t[i+1]) - known
        unew = levenberg(F,known)
        u[i+1] = unew[:,end]
    end
    return t,u
end

# This version is needed to work with scalar problems using levenberg().
function am2(dudt,tspan,u0::Number,n)
    f = (x,t) -> [dudt(x[1],t)]
    t,u = am2(f,tspan,[u0],n)
    u = [u[1] for u=u]
    return t,u
end

max_stars_repo_name: snowdj/fnc-extras

max_issues_repo_name: snowdj/fnc-extras

max_forks_repo_name: snowdj/fnc-extras


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_input(id="first", value = 1),
    dcc_input(id="second"),
    dcc_input(id="third")
end

callback!(
    ClientsideFunction("clientside", "add1_break_at_11"),
    app, Output("second", "value"), Input("first", "value")
    )

callback!(
    ClientsideFunction("clientside", "add1_break_at_11"),
    app, Output("third", "value"), Input("second","value")
    )

run_server(app)


text:
# FIXME: patches to improve inferrability of Base methods
include("patches.jl")

# test utilities are extracted into a separate file for easier interactive testing from REPL
# i.e. julia-dev -i test/interactive_utils.jl
include("interactive_utils.jl")

include(normpath(@__DIR__, "..", "benchmark", "benchmark_utils.jl"))

# stuff used across tests
using Test

const FIXTURE_DIR = normpath(@__DIR__, "fixtures")

const ERROR_REPORTS_FROM_SUM_OVER_STRING = let
    interp, frame = profile_call(sum, (String,))
    @test !isempty(interp.reports)
    interp.reports
end

function test_sum_over_string(ers)
    @test !isempty(ers)
    for target in ERROR_REPORTS_FROM_SUM_OVER_STRING
        @test any(ers) do er
            return er.msg == target.msg && er.sig == target.sig
        end
    end
end
test_sum_over_string(res::JET.VirtualProcessResult) =
    test_sum_over_string(res.inference_error_reports)
test_sum_over_string(interp::JETInterpreter) = test_sum_over_string(interp.reports)

@testset "JET.jl" begin
    @testset "virtualprocess.jl" begin
        include("test_virtualprocess.jl")
    end

    @testset "abstractinterpretation.jl" begin
        include("test_abstractinterpretation.jl")
    end

    @testset "tfuncs.jl" begin
        include("test_tfuncs.jl")
    end

    @testset "caching" begin
        include("test_caching.jl")
    end

    # tests with Windows-paths is just an hell
    @static Sys.iswindows() || @testset "print.jl" begin
        include("test_print.jl")
    end

    @testset "performance" begin
        include("test_performance.jl")
    end
end

# # test self profiling
# using JET
# interp, frame = let
#     m = first(methods(JET.profile_file))
#     interp = JET.JETInterpreter()
#     JET.profile_method!(interp, m)
#
#     # JET.print_reports(stdout::IO, interp.reports)
# end


text:
struct Nested{T}
    o::T
end

Base.iterate(u::Nested, i = u.o) = i isa AbstractObservable ? (i, i[]) : nothing

mutable struct Flatten <: AbstractObservable{Any}
    o::AbstractObservable
    list::Vector{Tuple{AbstractObservable, Function}}
    inner_id::String
    output::Observable{Any}
    pair::ObservablePair
    function Flatten(o)
        inner_obs = foldl((x,y)->y, Nested(o))
        inner_id = obsid(inner_obs)
        output = Observable{Any}(inner_obs[])
        p = ObservablePair(inner_obs, output)
        n = new(o, Tuple{AbstractObservable, Function}[], inner_id, output, p)
        for i in Nested(o)
            f = updater(i, n)
            on(f, i)
            push!(n.list, (i, f))
        end
        n
    end
end

flatten(x) = x
flatten(obs::AbstractObservable) = observe(Flatten(obs))

observe(o::Flatten) = o.output

function updater(i, f::Flatten)
    function (val)
        ind = findfirst(t -> obsid(t[1]) == obsid(i), f.list)
        if ind !== nothing
            for ii in length(f.list):-1:(ind+1)
                obs, func = pop!(f.list)
                off(obs, func)
            end
            for outer i in Nested(val)
                func = updater(i, f)
                on(func, i)
                push!(f.list, (i, func))
            end
            if obsid(i) != f.inner_id
                off(f.pair)
                f[] = i[]
                f.pair = ObservablePair(i, f.output)
                f.inner_id = obsid(i)
            end
        end
    end
end


text:
module JungleHelperGrablessGoldenBerry

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/TreeDepthController" GrablessGoldenBerry(x::Integer, y::Integer, nodes::Array{Tuple{Integer, Integer}, 1}=Tuple{Integer, Integer}[])

const placements = Ahorn.PlacementDict(
    "Golden Strawberry (Grabless) (Jungle Helper)" => Ahorn.EntityPlacement(
        GrablessGoldenBerry,
        "point"
    )
)

const sprite = "collectables/goldberry/wings01"
const seedSprite = "collectables/goldberry/seed00"

Ahorn.nodeLimits(entity::GrablessGoldenBerry) = 0, -1

function Ahorn.selection(entity::GrablessGoldenBerry)
    x, y = Ahorn.position(entity)

    nodes = get(entity.data, "nodes", ())

    res = Ahorn.Rectangle[Ahorn.getSpriteRectangle(sprite, x, y)]

    for node in nodes
        nx, ny = node

        push!(res, Ahorn.getSpriteRectangle(seedSprite, nx, ny))
    end

    return res
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::GrablessGoldenBerry)
    x, y = Ahorn.position(entity)

    for node in get(entity.data, "nodes", ())
        nx, ny = node

        Ahorn.drawLines(ctx, Tuple{Number, Number}[(x, y), (nx, ny)], Ahorn.colors.selection_selected_fc)
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::GrablessGoldenBerry, room::Maple.Room)
    x, y = Ahorn.position(entity)

    nodes = get(entity.data, "nodes", ())

    for node in nodes
        nx, ny = node

        Ahorn.drawSprite(ctx, seedSprite, nx, ny)
    end

    Ahorn.drawSprite(ctx, sprite, x, y)
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
include(joinpath(dirname(pathof(CImGui)), "..", "examples", "demo_window.jl"))

DoGui(() -> ShowDemoWindow(Ref(true)))


text:
"""
    flatten_levels(iter, lvls)

Flatten, or merge levels.

`flatten_levels(iter, 1)` coresponds to `Base.Iterators.flatten(iter)`

"""
function flatten_levels(iter, lvls)
    apply_at_level(Base.Iterators.flatten, iter, lvls)
end

"""
    consolidate_levels(iter, levels)

Applies `consolidate` on the given levels.
To turn those levels from iterators into something indexable.
`consolidate` is basically `collect` but it doesn't return `Vectors` just something indexable

See also `full_consolidate`, which is an alias for doing this on all levels.
"""
function consolidate_levels(iter, lvls)
    apply_at_level(consolidate, iter, lvls)
end

"""
    full_consolidate(iter)

Collect all iterators in the structor,
converting an Iterator of Iterators of ...
into a Vector of Vectors of ...
"""
full_consolidate(iter) =consolidate_levels(iter, ALL_LEVELS)

"""
    join_levels(iter, lvl_delims::Associative)

This is the generalisation of `join`.
Use this as `join_levels(doc_sents_words, Dict(1=>"\n", 2=" "))`,
to join each word with a space, and each sentence  line with a line.

Note that the elements in the level being joined are converted to strings via `string`.
This is find if they are already strings, or characters, or numbers etc.
But if they are not, e.g. they were nonmerged layers then this often going to be unusual output.
But what exactly did you expect it to do? (Answer? Throw an error)
"""
function join_levels(iter, lvl_delims::AbstractDict)
    apply_at_level(iter, Dict(lvl => (xs->join(xs,delim))  for (lvl,delim) in lvl_delims))
end


text:
# basic code for the 3D large-eddy simulations of a line plume

# imports - note that the Project.toml file needs to be present for these to work correctly
using DelimitedFiles, Printf, Interpolations
using Plots
using CuArrays, CUDAnative, CUDAdrv

using Oceananigans
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters
using Oceananigans.AbstractOperations
using Oceananigans.Utils
using Oceananigans: Face, Cell

#####
##### Preliminaries
#####

# some constants
const km = 1000
const Ω_Earth = 7.292115e-5  # [s⁻¹]
const φ = -75  # degrees latitude

# architecture
# modify as appropriate: if no GPU, set arch = CPU() and delete the following line
arch = GPU()
device!(CuDevice(2))

FT = Float64

# domain details (grid cells, then dimensions in m)
Nx = 512
Ny = 512
Nz = 96

Lx = 5km
Ly = 5km
Lz = 400

end_time = 6hour

g = 9.81
β = 7.8*10^(-4)

# SET UP INITIAL CONDITION
# Option 1: Make up salinity and temperature profiles based on prescribed N
#
#Nb = 0.003
#ΔS = (Lz*Nb^2)/(g*β)
#zC = collect(((-Lz:Lz/Nz:0).+Lz/(2*Nz))[1:end-1])
#T₀ = 1*ones(Nz)
#S₀ = reverse(collect(34.6-ΔS:ΔS/(Nz-1):34.6))

# Option 2: Read T and S data from file
zC = collect(((-Lz:Lz/Nz:0).+Lz/(2*Nz))[1:end-1])
zdata = (collect(1:1:5400))
T = (readdlm("pig_data/pig_outflow_2014_T.txt")[:])
S = (readdlm("pig_data/pig_outflow_2014_S.txt")[:])

# Get the indices of all the non-NaN values.
T_good_inds = findall(!isnan, T)
S_good_inds = findall(!isnan, S)

# Create T and S arrays that do not contain NaNs, along with corresponding
# z values.
T_good = T[T_good_inds]
S_good = S[S_good_inds]

z_T = zdata[T_good_inds]
z_S = zdata[S_good_inds]

# Linearly interpolate T and S profiles to model grid.
Ti = LinearInterpolation(z_T, T_good, extrapolation_bc=Interpolations.Flat())
Si = LinearInterpolation(z_S, S_good, extrapolation_bc=Interpolations.Flat())

T₀ = Ti.(-zC)
S₀ = Si.(-zC)

# end initial condition setup

# convert to CuArray if we are running on a GPU
if arch == GPU()
    T₀ = CuArray(T₀)
    S₀ = CuArray(S₀)
end

# Meltwater source location - implemented as a box
# width in terms of grid cells
widthx = 100
widthy = 10
source_corners = ((256-widthx/2,0,0),(256+widthx/2,widthy,1))

# specify the integrated buoyancy source (m^4/s^3) (i.e. F)
# then calculate salinity fluxes into bottom boxes
B_flux = 1
A = (Ly/Ny)*(source_corners[2][2] - source_corners[1][2])*(Lx/Nx)*(source_corners[2][1] - source_corners[1][1])
S_flux = -B_flux/(A*g*β*(Lz/Nz))

# meltwater passive tracer flux
m_flux = 1

# Forcing functions
@inline T_forcing(i, j, k, grid, time, U, C, p) =
    @inbounds ifelse(p.source_corners[1][1]<=i<=p.source_corners[2][1] && p.source_corners[1][2]<=j<=p.source_corners[2][2] && p.source_corners[1][3]<=k<=p.source_corners[2][3], 0, 0)

@inline S_forcing(i, j, k, grid, time, U, C, p) =
    @inbounds ifelse(p.source_corners[1][1]<=i<=p.source_corners[2][1] && p.source_corners[1][2]<=j<=p.source_corners[2][2] && p.source_corners[1][3]<=k<=p.source_corners[2][3], p.S_flux, 0)

@inline m_relax(i, j, k, grid, time, U, C, p) =
    @inbounds ifelse(p.source_corners[1][1]<=i<=p.source_corners[2][1] && p.source_corners[1][2]<=j<=p.source_corners[2][2] && p.source_corners[1][3]<=k<=p.source_corners[2][3], p.m_flux, 0)

params = (source_corners=source_corners, S_flux=S_flux, m_flux = m_flux, T₀=T₀,S₀=S₀)

forcing = ModelForcing(T = T_forcing, S = S_forcing, meltwater = m_relax)

#####
##### Set up model and simulation
#####

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

# equation of state - pick one!
# eos = LinearEquationOfState()
eos = RoquetIdealizedNonlinearEquationOfState(:freezing)

model = IncompressibleModel(
           architecture = arch,
             float_type = FT,
                   grid = grid,
                tracers = (:T, :S, :meltwater),
               coriolis = FPlane(rotation_rate=Ω_Earth, latitude=φ),
               buoyancy = SeawaterBuoyancy(equation_of_state=eos),
                closure = AnisotropicMinimumDissipation(),
                forcing = forcing,
             parameters = params
)

#####
##### Setting up initial conditions
#####

T₀_3D = repeat(reshape(T₀, 1, 1, Nz), Nx, Ny, 1)
S₀_3D = repeat(reshape(S₀, 1, 1, Nz), Nx, Ny, 1)

set!(model.tracers.T, T₀_3D)
set!(model.tracers.S, S₀_3D)

#####
##### Write out 3D fields and slices to NetCDF files.
#####

# Define vorticity computation
u, v, w = model.velocities
vorticity_operation = ∂x(v) - ∂y(u)
ω = Field(Face, Face, Cell, model.architecture, model.grid, TracerBoundaryConditions(grid))
vorticity_computation = Computation(vorticity_operation, ω)

function get_vorticity(model)
    compute!(vorticity_computation)
    return Array(interior(ω))
end

fields = Dict(
        "u" => model.velocities.u,
        "v" => model.velocities.v,
        "w" => model.velocities.w,
        "T" => model.tracers.T,
        "S" => model.tracers.S,
"meltwater" => model.tracers.meltwater,
       "nu" => model.diffusivities.νₑ,
   "kappaT" => model.diffusivities.κₑ.T,
   "kappaS" => model.diffusivities.κₑ.S,
   "vorticity" => get_vorticity
)

dimensions = Dict(
    "vorticity" => ("xF", "yF", "zC")
)

output_attributes = Dict(
"meltwater" => Dict("longname" => "Meltwater concentration"),
       "nu" => Dict("longname" => "Nonlinear LES viscosity", "units" => "m^2/s"),
   "kappaT" => Dict("longname" => "Nonlinear LES diffusivity for temperature", "units" => "m^2/s"),
   "kappaS" => Dict("longname" => "Nonlinear LES diffusivity for salinity", "units" => "m^2/s"),
   "vorticity" => Dict("longname" => "Vorticity", "units" => "1/s")
)

eos_name(::LinearEquationOfState) = "LinearEOS"
eos_name(::RoquetIdealizedNonlinearEquationOfState) = "RoquetEOS"
prefix = "$(Nx)x$(Ny)x$(Nz)_F$(B_flux)_w$(widthx)x$(widthy)c"

#####
##### Print banner
#####

@printf("""

    Simulating ocean dynamics of meltwater outflow from beneath Antarctic ice shelves
        N : %d, %d, %d
        L : %.3g, %.3g, %.3g [km]
        Δ : %.3g, %.3g, %.3g [m]
        φ : %.3g [latitude]
        f : %.3e [s⁻¹]
     days : %d
        F : %.2f [m⁴s⁻³]
  closure : %s
      EoS : %s

""", model.grid.Nx, model.grid.Ny, model.grid.Nz,
     model.grid.Lx / km, model.grid.Ly / km, model.grid.Lz / km,
     model.grid.Δx, model.grid.Δy, model.grid.Δz,
     φ, model.coriolis.f, end_time / day,
     B_flux,
     typeof(model.closure), typeof(model.buoyancy.equation_of_state))

#####
##### Time step!
#####

# Wizard utility that calculates safe adaptive time steps.
wizard = TimeStepWizard(cfl=0.2, Δt=1second, max_change=1.2, max_Δt=10second)

# Number of time steps to perform at a time before printing a progress
# statement and updating the adaptive time step.
Ni = 20

# CFL utilities for reporting stability criterions.
cfl = AdvectiveCFL(wizard)
dcfl = DiffusiveCFL(wizard)

function progress_statement(simulation)
    model = simulation.model
    C_mw = model.tracers.meltwater  # Convenient alias

    # Calculate simulation progress in %.
    progress = 100 * (model.clock.time / end_time)

    # Find maximum velocities.
    umax = maximum(abs, model.velocities.u.data.parent)
    vmax = maximum(abs, model.velocities.v.data.parent)
    wmax = maximum(abs, model.velocities.w.data.parent)

    # Find maximum ν and κ.
    νmax = maximum(model.diffusivities.νₑ.data.parent)
    κmax = maximum(model.diffusivities.κₑ.T.data.parent)

    # Print progress statement.
    i, t = model.clock.iteration, model.clock.time
    @printf("[%06.2f%%] i: %d, t: %5.2f days, umax: (%6.3g, %6.3g, %6.3g) m/s, CFL: %6.4g, νκmax: (%6.3g, %6.3g), νκCFL: %6.4g, next Δt: %8.5g s\n",
            progress, i, t / day, umax, vmax, wmax, cfl(model), νmax, κmax, dcfl(model), wizard.Δt)
end

_ω = get_vorticity(model)
#display(_ω)
#@show typeof(_ω)

# Simulation that manages time stepping.
simulation = Simulation(model, Δt=wizard, stop_time=end_time, progress=progress_statement, progress_frequency=Ni)

simulation.output_writers[:fields] =
    NetCDFOutputWriter(model, fields, filename = prefix * "fields.nc",
                       interval = hour, output_attributes = output_attributes,
		       dimensions = dimensions)

fields_without_vorticity = copy(fields)
delete!(fields_without_vorticity,"vorticity")

simulation.output_writers[:along_channel_slice] =
    NetCDFOutputWriter(model, fields_without_vorticity,
		       filename = prefix * "yz.nc",
                       interval = 5minute, output_attributes = output_attributes,
		       dimensions = dimensions, xC = Int(Nx/2), xF = Int(Nx/2))

simulation.output_writers[:along_front_slice] =
    NetCDFOutputWriter(model, fields_without_vorticity,
		       filename = prefix * "xz.nc",
                       interval = 5minute, output_attributes = output_attributes,
                       dimensions = dimensions, yC = 1, yF = 1)

run!(simulation)

for ow in simulation.output_writers
     ow isa NetCDFOutputWriter && close(ow)
end


text:
@testset "Input" begin
    #
    # check_input_fname
    #
    script1 = joinpath(J.PATHS[:assets], "script1.jl")
    write(script1, "1+1")
    fp, d, fn = J.check_input_rpath("script1.jl")
    @test fp == script1
    @test d == J.PATHS[:assets]
    @test fn == "script1"
    @test_throws ArgumentError J.check_input_rpath("script2.jl")

    #
    # resolve_lx_input_hlcode
    #
    r = J.resolve_lx_input_hlcode("script1.jl", "julia")
    r2 = J.resolve_lx_input_othercode("script1.jl", "julia")
    @test r == "
1+1
"
    @test r2 == r

    #
    # resolve_lx_input_plainoutput
    #
    mkpath(joinpath(J.PATHS[:assets], "output"))
    plain1 = joinpath(J.PATHS[:assets], "output", "script1.out")
    write(plain1, "2")

    r = J.resolve_lx_input_plainoutput("script1.jl")
    @test r == "
2
"
end

@testset "LX input" begin
    write(joinpath(J.PATHS[:assets], "s1.jl"), "println(1+1)")
    write(joinpath(J.PATHS[:assets], "output", "s1a.png"), "blah")
    write(joinpath(J.PATHS[:assets], "output", "s1.out"), "blih")
    st = raw"""
        Some string
        \input{julia}{s1.jl}
        Then maybe
        \input{output:plain}{s1.jl}
        Finally img:
        \input{plot:a}{s1.jl}
        done.
        """ * J.EOS;

    J.def_GLOBAL_PAGE_VARS!()
    J.def_GLOBAL_LXDEFS!()

    m, _ = J.convert_md(st, collect(values(J.GLOBAL_LXDEFS)))
    h = J.convert_html(m, J.PageVars())

    @test occursin("
Some string
$(read(joinpath(J.PATHS[:assets], "s1.jl"), String))
", h)
    @test occursin("Then maybe
$(read(joinpath(J.PATHS[:assets], "output", "s1.out"), String))
", h)
    @test occursin("Finally img:  done.", h)
end

@testset  "Input MD" begin
    mkpath(joinpath(J.PATHS[:assets], "ccc"))
    fp = joinpath(J.PATHS[:assets], "ccc", "asset1.md")
    write(fp, "blah **blih**")
    st = raw"""
        Some string
        \textinput{ccc/asset1}
        """ * J.EOS;
    @test isapproxstr(st |> conv, "
Some string blah blih
")
end


text:
using LinearAlgebra, Statistics
using Distributions, QuasiMonteCarlo, Random
using IterTools, ProgressBars
using LaTeXStrings
using PyPlot
using NPZ, JLD

SEED    = 42
Random.seed!(SEED)

# binary logistic regression objective with l2 regularization
a6a = npzread("a6a.npy")
y = a6a[:,1]
X = a6a[:,2:end]

RUNS    = 1
EPOCHS  = 1000
n       = size(X)[1]
d       = size(X)[2]
T       = n * EPOCHS
λ       = 1e-4

# Quasi-random distributions
qrngs = Dict([("IID Uniform"     , UniformSample()),
              ("Sobol"           , SobolSample())])
qrng_names = vcat(collect(keys(qrngs)), ["RR", "SO"])

function grad_errors(T::Int64, wtau::Vector{Float64},
                     fullgrad::Vector{Float64}, ordering::Vector{Int64})
    errs = zeros(T)
    acc = zeros(d)
    for m in 1:T
        i = ordering[m]
        xi = X[i,:]
        yi = y[i]
        res = - yi / (1 + exp(yi * xi'wtau))
        g = xi * res + λ * wtau
        acc +=  g
        errs[m] = norm(acc/m - fullgrad)^2
    end
    return errs
end

println("begin computing gradient errors...")
result = Dict()

w0 = zeros(d)
res = -y ./ (1 .+ exp.(y .* X * w0))
full_grad = reshape(sum(X .* res, dims = 1) ./ n, d) + λ * w0

for qrng_name in qrng_names
    err_list = zeros(RUNS, T)

    for run = 1:RUNS
        if qrng_name in collect(keys(qrngs))
            distrib = qrngs[qrng_name]
            # generate QMC indices
            sample_indices = [Int(ceil(i)) for i in QuasiMonteCarlo.sample(T, 0, n, distrib)]
        else
            if qrng_name == "SO"
                Random.seed!(SEED + run)
                perm = randperm!(collect(1:n))
                sample_indices = repeat(perm, EPOCHS)

            elseif qrng_name == "RR"
                sample_indices = zeros(Int64, T)
                for epoch = 1:EPOCHS
                    Random.seed!(SEED + run + epoch)
                    perm = randperm!(collect(1:n))
                    sample_indices[(epoch-1)*n+1:epoch*n] = perm
                end
            else
                ArgumentError("SGD variant does not exist.")
            end
        end

        err_list[run, :] = grad_errors(T, w0, full_grad, sample_indices)
    end
    result[qrng_name] = err_list
end

println("completed")
println("plotting...")

figure(figsize = (6,6))
colors = ["#1f77bf", "#ff7f0e", "#17becf", "#7f7f7f"]
for (q, qrng_name) in enumerate(["IID Uniform", "Sobol", "RR", "SO"])

    err_list = result[qrng_name]
    err_mean = mean(err_list, dims=1)'

    loglog(collect(1000:1000:T), err_mean[1000:1000:T], label=qrng_name, color=colors[q])
    ylim(bottom=1e-8)
end

legend(loc="upper right", fontsize=18)

xlabel(L"m", fontsize=20)
ylabel(L"\Vert \frac{1}{m}\Sigma_{t=0}^{m-1}\nabla R_n(w_0; x_t) - \nabla R_n(w_0)\Vert^2", fontsize=20)
title("Offline (a6a)", fontsize=24)
xticks(fontsize=20)
yticks(fontsize=20)

tight_layout()
plot_fn = joinpath(pwd(),"figs", "qmc-graderr-offline-a6a.pdf")
savefig(plot_fn)
println("plot saved at $plot_fn")

max_stars_repo_name: EugeneLYC/qmc-ordering

max_issues_repo_name: EugeneLYC/qmc-ordering

max_forks_repo_name: EugeneLYC/qmc-ordering


text:
### Transformations

abstract type AbstractDataTransform end

# apply the transform
"""
    transform!(t::AbstractDataTransform, x)

Apply transformation `t` to vector or matrix `x` in place.
"""
transform!(t::AbstractDataTransform, x::AbstractArray{<:Real,1}) = transform!(x, t, x)
transform!(t::AbstractDataTransform, x::AbstractArray{<:Real,2}) = transform!(x, t, x)

"""
    transform(t::AbstractDataTransform, x)

Return a row-standardized vector or matrix `x` using `t` transformation.
"""
transform(t::AbstractDataTransform, x::AbstractArray{<:Real,1}) = transform!(similar(x), t, x)
transform(t::AbstractDataTransform, x::AbstractArray{<:Real,2}) = transform!(similar(x), t, x)

# reconstruct the original data from transformed values
"""
    reconstruct!(t::AbstractDataTransform, y)

Perform an in-place reconstruction into an original data scale from a row-transformed
vector or matrix `y` using `t` transformation.
"""
reconstruct!(t::AbstractDataTransform, y::AbstractArray{<:Real,1}) = reconstruct!(y, t, y)
reconstruct!(t::AbstractDataTransform, y::AbstractArray{<:Real,2}) = reconstruct!(y, t, y)

"""
    reconstruct(t::AbstractDataTransform, y)

Return a reconstruction of an originally scaled data from a row-transformed vector
or matrix `y` using `t` transformation.
"""
reconstruct(t::AbstractDataTransform, y::AbstractArray{<:Real,1})  = reconstruct!(similar(y), t, y)
reconstruct(t::AbstractDataTransform, y::AbstractArray{<:Real,2})  = reconstruct!(similar(y), t, y)

"""
    Standardization (Z-score transformation)
"""
struct ZScoreTransform{T<:Real} <: AbstractDataTransform
    dim::Int
    mean::Vector{T}
    scale::Vector{T}

    function ZScoreTransform(d::Int, m::Vector{T}, s::Vector{T}) where T
        lenm = length(m)
        lens = length(s)
        lenm == d || lenm == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        lens == d || lens == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        new{T}(d, m, s)
    end
end

function Base.getproperty(t::ZScoreTransform, p::Symbol)
    if p === :indim || p === :outdim
        return t.dim
    else
        return getfield(t, p)
    end
end

"""
    fit(ZScoreTransform, X; center=true, scale=true)

Fit standardization parameters to `X` and return a `ZScoreTransform` transformation object.

# Arguments

* `data`: matrix  of samples to fit transformation parameters.

# Keyword arguments

* `center`: if `true` (the default) center data so that its mean is zero.

* `scale`: if `true` (the default) scale the data so that its variance is equal to one.

# Examples

```jldoctest
julia> using StatsBase

julia> X = [0.0 -0.5 0.5; 0.0 1.0 2.0]
2×3 Array{Float64,2}:
 0.0  -0.5  0.5
 0.0   1.0  2.0

julia> dt = fit(ZScoreTransform, X)
ZScoreTransform{Float64}(2, [0.0, 1.0], [0.5, 1.0])

julia> StatsBase.transform(dt, X)
2×3 Array{Float64,2}:
  0.0  -1.0  1.0
 -1.0   0.0  1.0
```
"""
function fit(::Type{ZScoreTransform}, X::AbstractArray{<:Real,2}; center::Bool=true, scale::Bool=true)
    d, n = size(X)
    n >= 2 || error("X must contain at least two columns.")

    T = eltype(X)
    m, s = mean_and_std(X, 2)

    return ZScoreTransform(d, (center ? vec(m) : zeros(T, 0)),
                              (scale ? vec(s) : zeros(T, 0)))
end

function transform!(y::AbstractVecOrMat{<:Real}, t::ZScoreTransform, x::AbstractVecOrMat{<:Real})
    d = t.dim
    size(x,1) == size(y,1) == d || throw(DimensionMismatch("Inconsistent dimensions."))
    n = size(y,2)
    size(x,2) == n || throw(DimensionMismatch("Inconsistent dimensions."))

    m = t.mean
    s = t.scale

    if isempty(m)
        if isempty(s)
            if x !== y
                copyto!(y, x)
            end
        else
            broadcast!(/, y, x, s)
        end
    else
        if isempty(s)
            broadcast!(-, y, x, m)
        else
            broadcast!((x,m,s)->(x-m)/s, y, x, m, s)
        end
    end
    return y
end

function reconstruct!(x::AbstractVecOrMat{<:Real}, t::ZScoreTransform, y::AbstractVecOrMat{<:Real})
    d = t.dim
    size(x,1) == size(y,1) == d || throw(DimensionMismatch("Inconsistent dimensions."))
    n = size(y,2)
    size(x,2) == n || throw(DimensionMismatch("Inconsistent dimensions."))

    m = t.mean
    s = t.scale

    if isempty(m)
        if isempty(s)
            if y !== x
                copyto!(x, y)
            end
        else
            broadcast!(*, x, y, s)
        end
    else
        if isempty(s)
            broadcast!(+, x, y, m)
        else
            broadcast!((y,m,s)->y*s+m, x, y, m, s)
        end
    end
    return x
end

"""
    Unit range normalization
"""
struct UnitRangeTransform{T<:Real}  <: AbstractDataTransform
    dim::Int
    unit::Bool
    min::Vector{T}
    scale::Vector{T}

    function UnitRangeTransform(d::Int, unit::Bool, min::Vector{T}, max::Vector{T}) where {T}
        lenmin = length(min)
        lenmax = length(max)
        lenmin == d || lenmin == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        lenmax == d || lenmax == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        new{T}(d, unit, min, max)
    end
end

function Base.getproperty(t::UnitRangeTransform, p::Symbol)
    if p === :indim || p === :outdim
        return t.dim
    else
        return getfield(t, p)
    end
end

# fit a unit transform
"""
    fit(UnitRangeTransform, X; center=true, scale=true)

Fit a scaling parameters to `X` and return transformation description.

# Arguments

* `data`: matrix  of samples to fit transformation parameters.

# Keyword arguments

* `center`: if `true` (the default) centere data around zero.

* `scale`: if `true` (the default) perform variance scaling.

# Examples

```jldoctest
julia> using StatsBase

julia> X = [0.0 -0.5 0.5; 0.0 1.0 2.0]
2×3 Array{Float64,2}:
 0.0  -0.5  0.5
 0.0   1.0  2.0

julia> dt = fit(UnitRangeTransform, X)
UnitRangeTransform{Float64}(2, true, [-0.5, 0.0], [1.0, 0.5])

julia> StatsBase.transform(dt, X)
2×3 Array{Float64,2}:
 0.5  0.0  1.0
 0.0  0.5  1.0
```
"""
function fit(::Type{UnitRangeTransform}, X::AbstractArray{<:Real,2}; unit::Bool=true)
    d, n = size(X)

    tmin = X[:, 1]
    tmax = X[:, 1]
    for j = 2:n
        @inbounds for i = 1:d
            if X[i, j] < tmin[i]
                tmin[i] = X[i, j]
            elseif X[i, j] > tmax[i]
                tmax[i] = X[i, j]
            end
        end
    end
    for i = 1:d
        @inbounds tmax[i] = 1 / (tmax[i] - tmin[i])
    end
    return UnitRangeTransform(d, unit, tmin, tmax)
end

function transform!(y::AbstractVecOrMat{<:Real}, t::UnitRangeTransform, x::AbstractVecOrMat{<:Real})
    d = t.dim
    size(x,1) == size(y,1) == d || throw(DimensionMismatch("Inconsistent dimensions."))
    n = size(x,2)
    size(y,2) == n || throw(DimensionMismatch("Inconsistent dimensions."))

    tmin = t.min
    tscale = t.scale

    if t.unit
        broadcast!((x,s,m) -> (x-m)*s, y, x, tscale, tmin)
    else
        broadcast!(*, y, x, tscale)
    end
    return y
end

function reconstruct!(x::AbstractVecOrMat{<:Real}, t::UnitRangeTransform, y::AbstractVecOrMat{<:Real})
    d = t.dim
    size(x,1) == size(y,1) == d || throw(DimensionMismatch("Inconsistent dimensions."))
    n = size(y,2)
    size(x,2) == n || throw(DimensionMismatch("Inconsistent dimensions."))

    tmin = t.min
    tscale = t.scale

    if t.unit
        broadcast!((y,s,m) -> y/s + m, x, y, tscale, tmin)
    else
        broadcast!(/, x, y, tscale)
    end
    return x
end

"""
    standardize(DT, X; kwargs...)

Return a row-standardized matrix `X` using `DT` transformation which is a subtype of `AbstractDataTransform`:

- `ZScoreTransform`
- `UnitRangeTransform`

# Example

```jldoctest
julia> using StatsBase

julia> standardize(ZScoreTransform, [0.0 -0.5 0.5; 0.0 1.0 2.0])
2×3 Array{Float64,2}:
  0.0  -1.0  1.0
 -1.0   0.0  1.0

julia> standardize(UnitRangeTransform, [0.0 -0.5 0.5; 0.0 1.0 2.0])
2×3 Array{Float64,2}:
 0.5  0.0  1.0
 0.0  0.5  1.0
```
"""
function standardize(::Type{DT}, X::AbstractArray{<:Real,2}; kwargs...) where {DT<:AbstractDataTransform}
    return transform(fit(DT, X; kwargs...), X)
end


text:
"""
This module allows users to push data using the [Pushbullet](http://pushbullet.com) service. It requires that the user have a Pushbullet account. The necessary Pushbullet token can be acquired from the [account settings](https://www.pushbullet.com/account) page.
"""
module Pushbullet

export user, devices, push_note, push_address, push_link, push_list

import HTTP
import JSON
using Compat

const PB_API_URL = "api.pushbullet.com/v2/"

struct PushbulletException <: Exception
    page :: AbstractString
    method :: Symbol
    status :: Int
end

Base.showerror(io :: IO, e :: PushbulletException) =
    print(io, uppercase(string(e.method)), " ", e.page, " HTTP Error Code: ", e.status)

function load_key(filepath="~/.pushbullet.key")
    fn = expanduser(filepath)
    open(fn, "r") do f
        chomp(readline(f))
    end
end

const PBKEY = load_key()

function api_call(page; method=:get, jsdata="")
    url = string("https://", PB_API_URL, page)

    header = @compat Dict("Access-Token" => PBKEY)

    if method == :get
        response = HTTP.get(url, header)
    elseif method == :post
        header["Content-Type"] = "application/json"
        response = HTTP.post(url, header, JSON.json(jsdata))
    else
        throw(ArgumentError("method must be :get or :post"))
    end

    if response.status == 200
        JSON.parse(String(response.body))
    else
        throw(PushbulletException(page, method, response.status))
    end
end

function user()
    api_call("users/me")
end

function matchattribute(device, key :: AbstractString, val :: Number)
    device[key] == val
end

function matchattribute(device, key :: AbstractString, val :: AbstractString)
    occursin(val, device[key])
end

function matchattribute(device, key :: AbstractString, val :: Regex)
    occursin(device[key], val)
end

function matchattribute(device, key :: AbstractString, predicate :: Function)
    predicate(device[key])
end

"""
Return a list of **active** devices on the supplied account which match the provided matchers. By default this returns a list of all active devices. Matchers may be provided as a *numeric value*, *string*, *regex*, or *predicate function* which will be compared ot the corresponding attribute of each device.

For example:

    Pushbullet.devices(nickname="nick")

This will return all devices with a nickname that is *exactly* `"nick"`.

We could instead call

    Pushbullet.devices(nickname=r"nick")

Which would return all devices whose nickname matches the regex `r"nick"`.

"""
function devices(;args...)
    devs = api_call("devices")["devices"]
    if isempty(args)
        filter!(dev -> dev["active"], devs)
    else
        filter!(devs) do dev
            dev["active"] && all(m -> matchattribute(dev, string(m[1]), m[2]), args)
        end
    end
end

"""
`iden` functions the same as `devices()` but returns on the `iden` attributes of the selected device objects.
"""
function iden(;args...)
    devs = devices(;args...)
    if isempty(devs)
        ""
    else
        devs[1]["iden"]
    end
end

function push(push_data)
    api_call("pushes", method=:post, jsdata=push_data)
end

function set_target!(push_data, target :: AbstractString)
    if isempty(target)
        push_data
    else
        push_data[:device_iden] = target
    end
end

function set_target!(push_data, target :: Dict{<:AbstractString, Any})
    set_target!(push_data, target["iden"])
end

function push_note(device_iden, title="", body="")
    push_data = @compat Dict(
        :type => "note",
        :title => title,
        :body => body)
    set_target!(push_data, device_iden)
    push(push_data)
end

function push_link(device_iden, title="", body="", url="")
    push_data = @compat Dict(
        :type => "link",
        :title => title,
        :body => body,
        :url => url)
    set_target!(push_data, device_iden)
    push(push_data)
end

function push_address(device_iden, name="", address="")
    push_data = @compat Dict(
        :type => "address",
        :name => name,
        :address => address)
    set_target!(push_data, device_iden)
    push(push_data)
end

function push_list(device_iden, title="", items=[])
    push_data = @compat Dict(
        :type => "list",
        :title => title,
        :items => items)
    set_target!(push_data, device_iden)
    push(push_data)
end

end # module


text:
using Test
using SHIPs, JuLIP, JuLIP.Testing, QuadGK, LinearAlgebra, SHIPs.JacobiPolys
using SHIPs: TransformedJacobi, transform, transform_d, eval_basis!,
             alloc_B, alloc_temp

##

@noinline get_IN(N) = collect((shpB.idx_Bll[N][1]+1):(shpB.idx_Bll[N][end]+shpB.len_Bll[N][end]))

# function barrier
@noinline gramian(N, shpB, Nsamples=100_000; normalise=false) =
   gramian(N, get_IN(N), shpB, Nsamples; normalise=normalise)

@noinline gramian(N, IN, shpB, Nsamples; normalise=false) =
   gramian(N, IN, alloc_temp(shpB), alloc_B(shpB), shpB, Nsamples, normalise)

function gramian(N, IN, tmp, B, shpB, Nsamples = 100_000, normalise = false)
   Zs = zeros(Int16, N)
   lenB = length(IN)
   G = zeros(Float64, lenB, lenB)
   for n = 1:Nsamples
      Rs = SHIPs.Utils.rand(shpB.J, N)
      eval_basis!(B, tmp, shpB, Rs, Zs, 0)
      for j = 1:lenB
         Bj = B[IN[j]]'
         @simd for i = 1:lenB
            @inbounds G[i,j] += Bj * B[IN[i]]
         end
      end
   end
   if normalise
      g = diag(G)
      for i = 1:lenB, j = 1:lenB
         G[i,j] /= sqrt(g[i]*g[j])
      end
   end
   return G / Nsamples
end

##
Nmax = 4
Nsamples = 1_000
rl, ru = 0.5, 3.0
fcut =  PolyCutoff2s(2, rl, ru)
trans = PolyTransform(2, 1.0)
shpB = SHIPBasis( SparseSHIP(6, 15), trans, fcut )

ctr = 0
nnz = 0
for N = 1:Nmax
   for νz in shpB.NuZ[N]
      kk, ll = SHIPs._kl(νz.ν, νz.izz, shpB.KL)
      for mm in SHIPs._mrange(ll)
         global ctr += 1
         global nnz += size( shpB.rotcoefs[N][ll], 2 )
      end
   end
end
length(shpB)
ctr
nnz

function fill_an_array!(A)
   for n = 1:length(A)
      A[n] = rand()
   end
   return A
end

A = zeros(100_000);
(@elapsed fill_an_array!(A)) * 1000

@show shpB.idx_Bll[end]
@show shpB.len_Bll[end]
@show length(shpB)

G = gramian(4, shpB, 10_000)
rank(G)
svdG = svd(G)
svdG.U

@info("Conditions numbers of gramians")
for N = 1:Nmax
   GN = gramian(N, shpB, Nsamples, normalise=false)
   @show N, cond(GN)
end

@info("Conditions numbers of normalised gramians")
for N = 1:Nmax
   GN = gramian(N, shpB, Nsamples, normalise=true)
   @show N, cond(GN)
end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
### push!

(s::Signal)((@nospecialize val)) = push!(s, val)
import Base.push!
function push!(s::Signal,(@nospecialize val), async::Bool = async_mode())
    try
        if s.strict_push
            strict_push!(s, val, async)
        else
            soft_push!(s, val, async)
        end
    catch e
        handle_err(e, stacktrace(catch_backtrace()))
    end
end

"""
    strict_push!(s::Signal, val, async = Signals.async_mode())

Set `s` to `val` and propagate into derived signals, is `async` is `true`(default)
then updates to derived signals will occur asynchronically.
"""
function strict_push!(s, (@nospecialize val), async = async_mode())
    if !async || isempty(pull_queue)
        soft_push!(s, val, async)
    else
        enqueue!(push_queue, (s,val))
        notify(eventloop_cond)
    end
    val
end
export strict_push!

function soft_push!(s,(@nospecialize val), async = async_mode())
    set_value!(s, val)
    propogate!(s, async)
    async && notify(eventloop_cond)
    val
end

function propogate!(s::Signal, async::Bool = async_mode())
    propagated(s) && return
    propagated(s, true)
    if isempty(s.children)
        pull_enqueue(s, async)
    else
        foreach(x -> propogate!(x, async), s.children)
    end
end

function pull_enqueue(s::Signal, async::Bool = async_mode())
    if async
        enqueue!(pull_queue, s)
    else
        pull!(s)
    end
end

### pull!

(s::Signal)() = pull!(s)

function pull!(s::Signal)
    if !valid(s)
        propagated(s, false)
        action(s)
    end
    return value(s)
end
pull!(@nospecialize x) = x

action(s::Signal) = s.action(s)

abstract type StandardPull <: PullType end
function (pa::PullAction{StandardPull, ARGS})(s::Signal) where ARGS
    args = pull_args(pa)
    if !valid(s)
        store!(s,pa.f(args...))
    end
    value(s)
end

nothing


text:
# usage example for package developer:
#
#     import ThresholdAPI: AbstractThresholdAlgorithm,
#                          find_threshold

"""
    AbstractThresholdAlgorithm <: AbstractImageFilter

The root type for `HistogramThreshold` package.

Any concrete threshold algorithm shall subtype it to support
[`find_threshold`](@ref) API.

# Examples

All algorithms in HistogramThreshold are called in the following pattern:

```julia
# Generate an algorithm instance.
f = Otsu()

# If you already have a histogram and concomitant edges at hand
# you can supply those directly.
t = find_threshold(histogram, edges, f)

# Alternatively, you can supply the raw data and specify histogram
# construction details.
t = find_threshold(data, f; nbins = 64)
```

For more examples, please check [`find_threshold`](@ref) and concret
algorithms.
"""
abstract type AbstractThresholdAlgorithm <: AbstractImageFilter end

function find_threshold(histogram::AbstractArray, edges::AbstractArray, f::AbstractThresholdAlgorithm)
    return f(histogram, edges)
end

function find_threshold(data::AbstractArray, f::AbstractThresholdAlgorithm ; nbins::Union{Int,Nothing} = 256)
    edges, counts  = build_histogram(data, nbins)
    #=
        The `counts` array stores at index 0 the frequencies that were below the
        first bin edge. Since we are seeking a threshold over the interval
        partitioned by `edges` we need to discard the first bin in `counts`
        so that the dimensions of `edges` and `counts` match.
    =#
    return f(view(counts, 1:lastindex(counts)), edges)
end

### Docstrings

"""
    find_threshold(data::AbstractArray, f::AbstractThresholdAlgorithm; nbins)
    find_threshold(histogram::AbstractArray, edges::AbstractArray, f::AbstractThresholdAlgorithm)

Find a suitable threshold in `data` using algorithm `f` upon constructing a histogram with `nbins`.
Instead of specifing the raw `data`, you can specify a histogram and accompanying edges directly.

# Output

A real number representing a threshold that can be used to split data into two parts.

# Examples

Just simply pass an algorithm to `find_threshold`:

```julia
using TestImages, HistogramThreshold
img = testimage("cameraman")
t = find_threshold(img, f ; nbins = 64)
```

```julia
using StatsBase, HistogramThreshold
data = vcat(ones(50,), zeros(50,))
h = fit(Histogram, data)
t = find_threshold(data, f ; nbins = 2)
```

"""
find_threshold

max_stars_repo_name: zygmuntszpak/HistogramThresholding

max_forks_repo_name: zygmuntszpak/HistogramThresholding