text:
# Builds a deterministic equivalent by joining the subproblems into
# a single MIP, then adding the constraints that there can only be
# a single investment along any path in the tree.
using DataStructures

struct DetEqModel
   problem::JuMP.Model
   tree::AbstractTree
   probabilities::Dict{AbstractTree,Float64}
   risk::Union{Risk,Array{Risk,1}}
end

"""
	DetEqModel(tree::AbstractTree,
               probabilities,
               sub_problem_builder::Function,
               solver
               discount_factor=1.0,
               risk=RiskNeutral,
               sideconstraints=nothing,
               parallel=false,
               check=true)

Define a deterministic equivalent model for the stochastic capacity expansion
problem.

### Required arguments
`tree` is a reference to a scenario tree

`probabilities` is either a function, which returns a dictionary of the probabilities
of all nodes in a tree, or simply the dictionary itself

`sub_problem_builder` is a function mapping a node to a JuMP model for each subproblems

`solver` is a reference to the optimizer used for this problem (with appropriate settings)

### Optional arguments
`discount_factor` is a number between 0 and 1 defining a constant discount factor along each arc
in the scenario tree

`risk` is a tuple with the two CVaR parameters: (λ, α)

`sideconstraints` is a function which specifies side constraints in the master problem, see
[Tutorial 9: Side-constraints](@ref) for further details.

`parallel` is a boolean, setting whether the sub-problems will be formulated in parallel

`check` is a boolean, which can be set to `false` to disable the validation of the JuDGE model.

### Examples
	deteq = DetEqModel(tree, ConditionallyUniformProbabilities, sub_problems,
                                    Gurobi.Optimizer)
	judge = DetEqModel(tree, probabilities, sub_problems, CPLEX.Optimizer,
                                    discount_factor=0.9, risk=(0.5,0.1)))
"""
function DetEqModel(tree::AbstractTree, probabilities, sub_problem_builder::Function, solver; discount_factor=1.0, risk::Union{Risk,Array{Risk,1}}=RiskNeutral(), sideconstraints=nothing, parallel=false, check=true)
   println("")
   println("Establishing deterministic equivalent model for tree: " * string(tree))
   if typeof(probabilities) <: Function
	   probabilities = probabilities(tree)
   end
   if typeof(probabilities)!=Dict{AbstractTree,Float64}
	   error("\'probabilities\' needs to be a dictionary mapping AbstractTree to Float64\nor a function that generates such a dictionary")
   end

   nodes=collect(tree)

   if parallel
	   sps=pmap(sub_problem_builder,nodes)
	   i=1
	   sub_problems=Dict()
	   for n in nodes
		   sub_problems[n]=sps[i]
		   i+=1
	   end
   else
	   sub_problems = Dict(i => sub_problem_builder(i) for i in nodes)
   end

   if check
	   print("Checking sub-problem format...")
	   check_specification_is_legal(sub_problems)
	   println("Passed")
   else
	   println("Skipping checks of sub-problem format")
   end
   if typeof(sub_problem_builder) <: Function
       JuDGE.scale_objectives(tree,sub_problems,discount_factor)
   end
   print("Building deterministic equivalent problem...")
   problem = build_deteq(sub_problems, tree, probabilities, solver, discount_factor, risk, sideconstraints)
   println("Complete")
   return DetEqModel(problem,tree,probabilities,risk)
end

function build_deteq(sub_problems::T where T <: Dict, tree::T where T <: AbstractTree, probabilities::Dict{AbstractTree,Float64}, solver, discount_factor::Float64,
	risk::Any, sideconstraints)

    model = JuMP.Model(solver)

    @objective(model,Min,0)

    model.ext[:vars]=Dict()
    model.ext[:master_vars]=Dict{AbstractTree,Dict{Symbol,Any}}()
    model.ext[:master_names]=Dict{AbstractTree,Dict{Symbol,Any}}()

    scen_con=Dict{Leaf,ConstraintRef}()
    scen_var=Dict{Leaf,VariableRef}()

	leafs=get_leafnodes(tree)
    for leaf in leafs
        scen_var[leaf]=@variable(model)
        scen_con[leaf]=@constraint(model,0==scen_var[leaf])
        set_objective_coefficient(model,scen_var[leaf],probabilities[leaf])
    end

	risk_objectives=AffExpr[]
	remain=1.0
	if typeof(risk)==Risk
		if (risk.α==1.0 || risk.λ==0.0) && risk.bound==nothing
			risk=[]
		else
			risk=[risk]
		end
	end
	for i in 1:length(risk)
		risk_objective=AffExpr(0.0)
		eta=@variable(model)
		for leaf in leafs
			v = @variable(model)
			w = @variable(model)
			set_lower_bound(v,0.0)
			set_lower_bound(w,0.0)
			if risk[i].offset==nothing
				@constraint(model,v>=eta-scen_var[leaf])
				@constraint(model,w>=scen_var[leaf]-eta)
				add_to_expression!(risk_objective,scen_var[leaf]*probabilities[leaf])
			else
				@constraint(model,v>=eta-scen_var[leaf]+risk[i].offset[leaf])
				@constraint(model,w>=scen_var[leaf]-risk[i].offset[leaf]-eta)
				add_to_expression!(risk_objective,(scen_var[leaf]-risk[i].offset[leaf])*probabilities[leaf])
			end
			add_to_expression!(risk_objective,probabilities[leaf]*(v+w/risk[i].α*(1-risk[i].α)))
		end
		remain-=risk[i].λ
		push!(risk_objectives,risk_objective)
	end

    for (node,sp) in sub_problems
        leafnodes=get_leafnodes(node)
        model.ext[:vars][node] = Dict()
        for variable in all_variables(sp)
            model.ext[:vars][node][variable] = JuDGE.copy_variable!(model, variable)
            if variable==sp.ext[:objective]
                for leaf in leafnodes
                    set_normalized_coefficient(scen_con[leaf],model.ext[:vars][node][variable],1.0)
                end
            end
        end
		copy_values=false
        loct=list_of_constraint_types(sp)

        for ct in loct
            for con in all_constraints(sp,ct[1],ct[2])
                con_obj=JuMP.constraint_object(con)
                if typeof(con_obj.func)==VariableRef
                    LHS=AffExpr(0.0)
                    add_to_expression!(LHS,1,model.ext[:vars][node][JuMP.constraint_object(con).func])
                elseif typeof(con_obj.func) <: GenericAffExpr
                    LHS=AffExpr(0.0)
                    for (v,c) in con_obj.func.terms
                        add_to_expression!(LHS,c,model.ext[:vars][node][v])
                    end
                elseif typeof(con_obj.func) <: GenericQuadExpr
                    LHS=QuadExpr(AffExpr(0.0),OrderedDict{UnorderedPair{VariableRef},Float64}())
                    for (v,c) in con_obj.func.terms
                        LHS.terms[UnorderedPair{VariableRef}(model.ext[:vars][node][v.a],model.ext[:vars][node][v.b])]=c
                    end
                    for (v,c) in con_obj.func.aff.terms
                        add_to_expression!(LHS.aff,c,model.ext[:vars][node][v])
                    end
                elseif typeof(con_obj.func) == Array{GenericAffExpr{Float64,VariableRef},1}
                    group=Array{AffExpr,1}()
                    for aff in con_obj.func
                        LHS=AffExpr(0.0)
                        for (v,c) in aff.terms
                            add_to_expression!(LHS,c,model.ext[:vars][node][v])
                        end
                        push!(group,LHS)
                    end
                else
                    error("Unsupported constraint type found: "*string(typeof(con_obj.func)))
                end
                set=con_obj.set
                if typeof(set)==MathOptInterface.GreaterThan{Float64}
                    @constraint(model,LHS>=set.lower)
                elseif typeof(set)==MathOptInterface.LessThan{Float64}
                    @constraint(model,LHS<=set.upper)
                elseif typeof(set)==MathOptInterface.EqualTo{Float64}
                    @constraint(model,LHS==set.value)
                elseif typeof(set)==MathOptInterface.SecondOrderCone
                    @constraint(model,group in SecondOrderCone())
                elseif typeof(set)==MathOptInterface.IndicatorSet{MathOptInterface.ACTIVATE_ON_ZERO,MathOptInterface.EqualTo{Float64}}
                    @constraint(model,!collect(keys(group[1].terms))[1] => {group[2]==set.set.value})
                elseif typeof(set)==MathOptInterface.IndicatorSet{MathOptInterface.ACTIVATE_ON_ZERO,MathOptInterface.LessThan{Float64}}
                    @constraint(model,!collect(keys(group[1].terms))[1] => {group[2]<=set.set.value})
                elseif typeof(set)==MathOptInterface.IndicatorSet{MathOptInterface.ACTIVATE_ON_ZERO,MathOptInterface.GreaterThan{Float64}}
                    @constraint(model,!collect(keys(group[1].terms))[1] => {group[2]>=set.set.value})
                elseif typeof(set)==MathOptInterface.IndicatorSet{MathOptInterface.ACTIVATE_ON_ONE,MathOptInterface.EqualTo{Float64}}
                    @constraint(model,collect(keys(group[1].terms))[1] => {group[2]==set.set.value})
                elseif typeof(set)==MathOptInterface.IndicatorSet{MathOptInterface.ACTIVATE_ON_ONE,MathOptInterface.LessThan{Float64}}
                    @constraint(model,collect(keys(group[1].terms))[1] => {group[2]<=set.set.value})
                elseif typeof(set)==MathOptInterface.IndicatorSet{MathOptInterface.ACTIVATE_ON_ONE,MathOptInterface.GreaterThan{Float64}}
                    @constraint(model,collect(keys(group[1].terms))[1] => {group[2]>=set.set.value})
                elseif typeof(set)!=MathOptInterface.ZeroOne && typeof(set)!=MathOptInterface.Integer
                    error("Unsupported constraint type found: "*string(typeof(set)))
                else
                    continue
                end
            end
        end
    end

    for (node,sp) in sub_problems
        model.ext[:master_vars][node] = Dict{Symbol,Any}()
        model.ext[:master_names][node] = Dict{Symbol,Any}()
        for (name,exps) in sp.ext[:expansions]
            if isa(exps,VariableRef)
                variable=sp.ext[:expansions][name]
                model.ext[:master_vars][node][name] =JuDGE.copy_variable!(model, variable)
                model.ext[:master_names][node][name] = string(variable)
            elseif typeof(exps) <: AbstractArray
                variables=sp.ext[:expansions][name]
                model.ext[:master_vars][node][name]=Dict()
                model.ext[:master_names][node][name]=Dict()
                for index in keys(exps)
					key=densekey_to_tuple(index)
                    model.ext[:master_vars][node][name][key] =JuDGE.copy_variable!(model, variables[index])
                    model.ext[:master_names][node][name][key]=string(variables[index])
                end
            end
			if sp.ext[:options][name][5]!=nothing
				if typeof(exps) <: AbstractArray
					for index in keys(exps)
						key=densekey_to_tuple(index)
						set_lower_bound(model.ext[:master_vars][node][name][key],sp.ext[:options][name][5])
					end
				else
					set_lower_bound(model.ext[:master_vars][node][name],sp.ext[:options][name][5])
				end
			end
			if sp.ext[:options][name][6]!=nothing
				if typeof(variable) <: AbstractArray
					for index in keys(exps)
						key=densekey_to_tuple(index)
						set_upper_bound(model.ext[:master_vars][node][name][key],sp.ext[:options][name][6])
					end
				else
					set_upper_bound(model.ext[:master_vars][node][name],sp.ext[:options][name][6])
				end
			end
        end
    end

    for leaf in get_leafnodes(tree)
        nodes=history(leaf)
        for n in eachindex(nodes)
            node=nodes[n]
            sp=sub_problems[node]
            df=discount_factor^depth(node)
            for (name,exps) in sp.ext[:expansions]
                interval=max(1,n-sp.ext[:options][name][3]-sp.ext[:options][name][2]+1):n-sp.ext[:options][name][2]
				disc=Dict{Int64,Float64}()
				for i in interval
					disc[i]=df/discount_factor^(i-1)
				end
                if isa(exps,VariableRef)
                    variable=sp.ext[:expansions][name]
                    cost_coef=df*coef(sp.ext[:capitalcosts],variable)
					for j in interval
						cost_coef+=disc[j]*coef(sp.ext[:ongoingcosts],variable)
					end
                    set_normalized_coefficient(scen_con[leaf],model.ext[:master_vars][node][name],cost_coef)
                elseif typeof(exps) <: AbstractArray
                    variables=sp.ext[:expansions][name]
                    for index in keys(exps)
						key=densekey_to_tuple(index)
                        cost_coef=df*coef(sp.ext[:capitalcosts],variables[index])
						for j in interval
							cost_coef+=disc[j]*coef(sp.ext[:ongoingcosts],variables[index])
						end
                        set_normalized_coefficient(scen_con[leaf],model.ext[:master_vars][node][name][key],cost_coef)
                    end
                end
            end
        end
    end

    for (node,sp) in sub_problems
        past=history(node)
        for (name,exps) in sp.ext[:expansions]
            interval=sp.ext[:options][name][2]+1:min(sp.ext[:options][name][2]+sp.ext[:options][name][3],length(past))
            if isa(exps,VariableRef)
                if sp.ext[:options][name][1]==:shutdown
                    @constraint(model,model.ext[:vars][node][exps]>=sum(model.ext[:master_vars][past[index]][name] for index in interval))
                elseif sp.ext[:options][name][1]==:expansion
                    @constraint(model,model.ext[:vars][node][exps]<=sum(model.ext[:master_vars][past[index]][name] for index in interval))
				elseif sp.ext[:options][name][1]==:enforced
					@constraint(model,model.ext[:vars][node][exps]==sum(model.ext[:master_vars][past[index]][name] for index in interval))
				elseif sp.ext[:options][name][1]==:state
					if node.parent==nothing
						@constraint(model, model.ext[:vars][node][exps] == -sp.ext[:options][name][7] + model.ext[:master_vars][node][name])
					else
						@constraint(model, model.ext[:vars][node][exps] == model.ext[:master_vars][node][name] - model.ext[:master_vars][node.parent][name])
					end
                end
                # if typeof(node)==Leaf && sp.ext[:options][name][1]==:shutdown
                #     @constraint(model,sum(model.ext[:master_vars][n][name] for n in history_function(node))<=1)
                # end
            elseif typeof(exps) <: AbstractArray
                for i in keys(exps)
					key=densekey_to_tuple(i)
                    if sp.ext[:options][name][1]==:shutdown
                        @constraint(model,model.ext[:vars][node][exps[i]]>=sum(model.ext[:master_vars][past[index]][name][key] for index in interval))
                    elseif sp.ext[:options][name][1]==:expansion
                        @constraint(model,model.ext[:vars][node][exps[i]]<=sum(model.ext[:master_vars][past[index]][name][key] for index in interval))
					elseif sp.ext[:options][name][1]==:enforced
                        @constraint(model,model.ext[:vars][node][exps[i]]==sum(model.ext[:master_vars][past[index]][name][key] for index in interval))
					elseif sp.ext[:options][name][1]==:state
						if node.parent==nothing
							@constraint(model, model.ext[:vars][node][exps[i]] == -sp.ext[:options][name][7] + model.ext[:master_vars][node][name][i])
						else
							@constraint(model, model.ext[:vars][node][exps[i]] == model.ext[:master_vars][node][name][i] - model.ext[:master_vars][node.parent][name][i])
						end
                    end
                    # if typeof(node)==Leaf && sp.ext[:options][name][1]==:shutdown
                    #     @constraint(model,sum(model.ext[:master_vars][n][name][i] for n in history_function(node))<=1)
                    # end
                end
            end
        end
    end

    if typeof(sideconstraints) <: Function
        map(Main.eval,unpack_expansions(model.ext[:master_vars])) #bring expansion variables into global scope
        sideconstraints(model,tree)
        map(Main.eval,clear_expansions(model.ext[:master_vars]))
    end

	objective_fn=objective_function(model)*remain

	for i in 1:length(risk_objectives)
		objective_fn+=risk_objectives[i]*risk[i].λ
		if risk[i].bound!=nothing
			surplus=@variable(model)
			set_lower_bound(surplus,0)
			@constraint(model,risk_objectives[i]<=risk[i].bound+surplus)
			objective_fn+=risk[i].penalty*surplus
		end
	end

	set_objective_function(model,objective_fn)

    model.ext[:scenario_obj]=scen_var

    return model
end

"""
	solve(deteq::DetEqModel)

Solve a determinisitc equivalent model.

### Required Arguments
`deteq` is the determinisitc equivalent model that we wish to solve.

### Example
    JuDGE.solve(deteq)
"""
function solve(deteq::DetEqModel)
    print("Solving deterministic equivalent formulation...")
    optimize!(deteq.problem)
    if termination_status(deteq.problem) == MathOptInterface.OPTIMAL
        println("Solved.")
    else
        println("Not solved: " * string(termination_status(deteq.problem)))
    end
end

function get_objval(deteq::DetEqModel; risk=deteq.risk)
	scenario_objs=Dict{Leaf,Float64}()

	for (leaf, var) in deteq.problem.ext[:scenario_obj]
		scenario_objs[leaf]=JuMP.value(var)
	end

	compute_objval(scenario_objs, deteq.probabilities, risk)
end


text:
#!/bin/bash
#=
export JULIA_PROJECT=$(dirname $(dirname $(readlink -e "${BASH_SOURCE[0]}")))
exec julia --color=yes --startup-file=no "${BASH_SOURCE[0]}" "$@"
=#

using HDF5

function getantnames(bfr)
    h5open(bfr) do h5
        h5["telinfo/antenna_names"][]
    end
end

function mkobsinfo(io::IO, antennas, instrument="BLUSE")
    println(io, """
# Name of instrument
instrument: BLUSE

# Input map.  This lists which antenna-polarizations were connected to which
# correlator inputs.  For dual-polarization inputs, the two polarization for a
# given input are flattened into two sequential elements, with the "first"
# polarization of the input given before the "second" polarization.  The
# antenna-polarization for an input is given as a two element array.  The first
# element specifies the antenna.  This can be given as a string for antenna
# name or an integer for antenna number.  The second element specifies the
# polarization.  Valid values are L,R,X,Y (case-insensitive).  Polarizations
# must be specified for each input, even for single polarization instruments.
input_map: ["""
    )

    for a in antennas
        println(io, "  [$a, x], [$a, y],")
    end

    println(io, "]")
end

if abspath(PROGRAM_FILE) == @__FILE__
    for arg in ARGS
        if isdir(arg)
            # Assume/require all files in a subdir to have the same obsinfo
            bfrnames = filter(endswith(".bfr5"), readdir(arg, join=true))
            if isempty(bfrnames)
                @warn "$arg has no bfr5 files" _module=nothing _file=nothing
                continue
            end
            bfrname = bfrnames[1]
        else
            if !ispath(arg)
                @warn "$arg does not exist" _module=nothing _file=nothing
                continue
            end
            bfrname = arg
        end

        obsinfoname = joinpath(dirname(bfrname), "obsinfo.yml")
        @info "$bfrname => $obsinfoname"
        open(obsinfoname, "w") do obsinfo
            mkobsinfo(obsinfo, getantnames(bfrname))
        end
    end
end

ext: jl


text:
using Test
using MeshCat
using FileIO
using CoordinateTransformations

# Demonstrate a couple of different lighting configurations, and
# show how the lighting settings can be manipulated programmatically.
@testset "Lighting Scenes" begin
    cat_mesh_path = joinpath(@__DIR__, "data", "meshes", "cat.obj")
    cat = load(cat_mesh_path)

    function add_cats!(vis::Visualizer)
        setobject!(vis["floor"], Rect(Vec(-10, -10, -0.01), Vec(20, 20, 0.01)))
        for x in range(-5, stop=5, step=2)
            for y in range(-5, stop=5, step=2)
                setobject!(vis["cat"][string(x)][string(y)], cat)
                settransform!(vis["cat"][string(x)][string(y)], Translation(x, y, 0) ∘ LinearMap(RotZ(π)) ∘ LinearMap(RotX(π/2)))
            end
        end
    end

    @testset "Default" begin
        vis = Visualizer()
        if !haskey(ENV, "CI")
            open(vis)
            wait(vis)
        end
        add_cats!(vis)
    end

    @testset "SpotLight with shadows" begin
        vis = Visualizer()
        if !haskey(ENV, "CI")
            open(vis)
            wait(vis)
        end
        add_cats!(vis)
        setprop!(vis["/Lights/SpotLight"], "visible", true)
        # To understand why we need the <object> here, see the documentation for
        # set_property at https://github.com/rdeits/meshcat
        setprop!(vis["/Lights/SpotLight/<object>"], "castShadow", true)
        setprop!(vis["/Lights/DirectionalLight"], "visible", false)
        # To understand why we need the <object> here, see the documentation for
        # set_property at https://github.com/rdeits/meshcat
        setprop!(vis["/Lights/AmbientLight/<object>"], "intensity", 0.44)
    end
end

@testset "Background" begin
    @testset "Background color" begin
        vis = Visualizer()
        if !haskey(ENV, "CI")
            open(vis)
            wait(vis)
        end
        setprop!(vis["/Background"], "top_color", colorant"red")
    end

    @testset "Background hidden" begin
        vis = Visualizer()
        if !haskey(ENV, "CI")
            open(vis)
            wait(vis)
        end
        setprop!(vis["/Background"], "visible", false)
    end

end


text:
# See p. 280 of Murphy's Machine Learning
# x_k1 = x_k - alpha * gr + mu * (x - x_previous)

immutable MomentumGradientDescent{L<:Function} <: Optimizer
    mu::Float64
    linesearch!::L
end

#= uncomment for v0.8.0
MomentumGradientDescent(; mu::Real = 0.01, linesearch = LineSearches.hagerzhang!) =
  MomentumGradientDescent(Float64(mu), linesearch)
=#

function MomentumGradientDescent(; mu::Real = 0.01, linesearch! = nothing,
                                   linesearch = LineSearches.hagerzhang!)
    linesearch = get_linesearch(linesearch!, linesearch)
    MomentumGradientDescent(Float64(mu), linesearch)
end

type MomentumGradientDescentState{T}
    @add_generic_fields()
    x_previous::Array{T}
    g::Array{T}
    f_x_previous::T
    s::Array{T}
    @add_linesearch_fields()
end

function initial_state{T}(method::MomentumGradientDescent, options, d, initial_x::Array{T})
    g = similar(initial_x)
    f_x = d.fg!(initial_x, g)

    MomentumGradientDescentState("Momentum Gradient Descent",
                         length(initial_x),
                         copy(initial_x), # Maintain current state in state.x
                         f_x, # Store current f in state.f_x
                         1, # Track f calls in state.f_calls
                         1, # Track g calls in state.g_calls
                         0, # Track h calls in state.h_calls
                         copy(initial_x), # Maintain current state in state.x_previous
                         g, # Store current gradient in state.g
                         T(NaN), # Store previous f in state.f_x_previous
                         similar(initial_x), # Maintain current search direction in state.s
                         @initial_linesearch()...) # Maintain a cache for line search results in state.lsr
end

function update_state!{T}(d, state::MomentumGradientDescentState{T}, method::MomentumGradientDescent)
    lssuccess = true
    # Search direction is always the negative gradient
    @simd for i in 1:state.n
        @inbounds state.s[i] = -state.g[i]
    end

    # Refresh the line search cache
    dphi0 = vecdot(state.g, state.s)
    LineSearches.clear!(state.lsr)
    push!(state.lsr, zero(T), state.f_x, dphi0)

    # Determine the distance of movement along the search line
    try
        state.alpha, f_update, g_update =
        method.linesearch!(d, state.x, state.s, state.x_ls, state.g_ls, state.lsr,
                           state.alpha, state.mayterminate)
        state.f_calls, state.g_calls = state.f_calls + f_update, state.g_calls + g_update
    catch ex
        if isa(ex, LineSearches.LineSearchException)
            lssuccess = false
            state.f_calls, state.g_calls = state.f_calls + ex.f_update, state.g_calls + ex.g_update
            state.alpha = ex.alpha
            Base.warn("Linesearch failed, using alpha = $(state.alpha) and exiting optimization.")
        else
            rethrow(ex)
        end
    end

    # Update current position
    @simd for i in 1:state.n
        # Need to move x into x_previous while using x_previous and creating "x_new"
        @inbounds tmp = state.x_previous[i]
        @inbounds state.x_previous[i] = state.x[i]
        @inbounds state.x[i] = state.x[i] + state.alpha * state.s[i] + method.mu * (state.x[i] - tmp)
    end
    (lssuccess == false) # break on linesearch error
end


text:
function refl_DT(t,β,ρ)

	nmed = 1.35
	ndet = 1.45
	μa = β[1]
    μsp = β[2]
    n = 1.35
    c = 29.9792458 # Speed of light cm/ns
    v = c/n # Speed of light in medium
	D = 1/(3*(μa + μsp)) # Diffusion Coefficient
	z0 = 1/μsp

	if (nmed == ndet)
		Afac = 1
	elseif nmed > ndet
		costhetac = sqrt(1 - (ndet/nmed).^2)
		R0 = ((nmed/ndet-1)/(nmed/ndet +1)).^2
		Afac = (2/(1-R0) -1 + abs(costhetac.^3))./(1-abs(costhetac.^2))
	else #nmed < ndet
		R0 = ((nmed/ndet-1)/(nmed/ndet +1)).^2
		Afac = 2/(1-R0) - 1
	end

	ze = 2*Afac*D
	zp = z0 + 2*ze

	Rt1 = @. -exp(-(ρ^2)/(4*D*v*t)-μa*v*t)/(2*t^(5/2)*(4*π*D*v)^3/2)
	Rt2 = @. (z0*exp(-(z0^2)/(4*D*v*t)) - zp*exp(-(zp^2)/(4*D*v*t)))
	Rt = @. Rt1*Rt2
    #Rt = @. (v/(4π*D*v*t)^1.5)*exp(-(ρ^2)/(4D*v*t)-μa*v*t)

    replace!(Rt, NaN => 0)
    m = findmax(Rt)
    Rt = Rt./m[1]
end

function refl_DT1(t,β,ρ)

	nmed = 1.35
	ndet = 1.45
	μa = β[1]
	μsp = β[2]
	#ρ = β[3]

	if (nmed == ndet)
		Afac = 1
	elseif nmed > ndet
		costhetac = sqrt(1 - (ndet/nmed).^2)
		R0 = ((nmed/ndet-1)/(nmed/ndet +1)).^2
		Afac = (2/(1-R0) -1 + abs(costhetac.^3))./(1-abs(costhetac.^2))
	else nmed < ndet
		R0 = ((nmed/ndet-1)/(nmed/ndet +1)).^2
		Afac = 2/(1-R0) - 1
	end

	z0 = 1/(μa+μsp)
	D = 1/(3*(μa + μsp))
	zb = 2*Afac*D

	v = 29.9792345/nmed # speed of light in medium

	Rt1 = @. v*exp(-μa*v*t)
	Rt1 = @. Rt1/((4*π*D*v*t)^1.5)
	Rt1 = @. Rt1*(exp(-(z0^2 + ρ^2)/(4*D*v*t)) - exp(-((2*zb + z0)^2 + ρ^2)/(4*D*v*t)))

	Rt2 = @. 3*exp(-μa*v*t)/(2*((4*π*D*v)^1.5)*(t^2.5))
	Rt2 = @. Rt2*(z0*exp(-((z0^2 + ρ^2)/(4*D*v*t))) + (2*zb + z0)*exp(-((2*zb+z0)^2 + ρ^2)/(4*D*v*t)))

	Rt = @. abs(Rt1)+abs(Rt2)
	replace!(Rt, NaN => 0)
	m = findmax(Rt)
	Rt = Rt./m[1]
end

max_stars_repo_name: heltonmc/load-asc-data

max_issues_repo_name: heltonmc/load-asc-data

max_forks_repo_name: heltonmc/load-asc-data


text:
###
### (::Tensor)(::TVector, :, TCovector, ...)
###

function (t::Tensor{<:Any,N})(v::Union{TVector,TCovector}; dims=N) where {N}
    N == 0 && throw(ArgumentError("Cannot contract Scalar with $(typeof(v))"))
    position(index_pos(t), dims) !== position(index_pos(v), 1) ||
        throw(ArgumentError("Tensor indices don't match"))

    # TODO: use an external package for tensor reductions
    parent = mapreduce(
        (x, y) -> x .* y, (x, y) -> x .+ y, eachslice(t.parent, dims=dims), v;
        init=zeros(promote_type(eltype(t), eltype(v)), size(t)[setdiff(1:end, dims)]),
    )
    deleteat(index_pos(t), dims)(parent)
end
(t::Tensor)(::Colon; dims=0) = t
function (t::Tensor{<:Any,N})(args...) where {N}
    foldl((t, (i,v)) -> t(v, dims=N-i+1), enumerate(args), init=t)
end

###
### raise/lower
###

function raise(t::Tensor{<:Any,N,ipos}, i::Int) where {N,ipos}
    IndexPos{N}(ipos.x | (UInt(1) << (i-1)))(t.parent)
end
raise(t::Tensor, inds::Vararg{Int}) = foldl((t, i) -> raise(t, i), inds, init=t)

function lower(t::Tensor{<:Any,N,ipos}, i::Int) where {N,ipos}
    IndexPos{N}(ipos.x & ~(UInt(1) << (i-1)))(t.parent)
end
lower(t::Tensor, inds::Vararg{Int}) = foldl((t, i) -> lower(t, i), inds, init=t)


text:
# Solve a SAA of the resource allocation model given scenarios of the demands

function solveSAAModel(demand_scen::Array{Float64,2})

	numScenarios::Int64 = size(demand_scen,1)

	mod = Model(solver=GurobiSolver(gurobi_env,Presolve=0,OutputFlag=0))

	@variable(mod, z[1:numResources] >= 0)
	@variable(mod, v[1:numResources,1:numCustomers,1:numScenarios] >= 0)
	@variable(mod, w[1:numCustomers,1:numScenarios] >= 0)

	@objective(mod, Min, sum(costVector[i]*z[i] for i=1:numResources) + (1.0/numScenarios)*sum(recourseCostCoeff[j]*w[j,s] for j=1:numCustomers for s=1:numScenarios))

	@constraint(mod, [i=1:numResources, s=1:numScenarios], sum(v[i,j,s] for j = 1:numCustomers) <= rho[i]*z[i])
	@constraint(mod, [j=1:numCustomers, s=1:numScenarios], sum(mu[i,j]*v[i,j,s] for i = 1:numResources) + w[j,s] >= demand_scen[s,j])

	status = solve(mod)
	z_soln = getvalue(z)
	objValue_soln::Float64 = getobjectivevalue(mod)

	return z_soln, objValue_soln
end

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_name: rohitkannan/DD-SAA


text:
using JLD2, Test

fn = joinpath(mktempdir(), "test.jld")
f = jldopen(fn, "w")
write(f, "/test_group_1/x1", 1)
@test !isempty(read(f, "test_group_1"))
write(f, "test_group_1/x2", 2)
write(f, "/test_group_2/x1", 3)
write(f, "test_group_3/x1", 4)
read(f, "test_group_3")["contained_group/x1"] = 5
JLD2.Group(f, "test_group_4")["x1"] = 6
JLD2.Group(read(f, "test_group_4"), "contained_group")["x1"] = 7
@test isempty(JLD2.Group(f, "empty_group"))
JLD2.Group(f, "group_with_one_group_child/empty_group")
@test !isempty(read(f, "group_with_one_group_child"))
@test_throws ArgumentError JLD2.Group(f, "test_group_4")
@test_throws ArgumentError write(f, "test_group_4/contained_group", 0)
@test !haskey(f, "path_that_does_not_exist")
@test !haskey(f, "/test_group_1/x1/path/that/does/not/exist")
@test !haskey(f["test_group_1"], "/test_group_1/x1/path/that/does/not/exist")
@test !haskey(f["test_group_1"], "x1/path/that/does/not/exist")
@test haskey(f, "/test_group_1/x1")
close(f)
f = jldopen(fn, "a")
@test haskey(f, "test_group_1")
@test haskey(f, "test_group_1/x1")

@test get(f, "path_that_does_not_exist", 2.0) == 2.0
@test get!(f, "path_created_by_get!", "x") == "x"
@test f["path_created_by_get!"] == "x"
@test get(JLD2.Group(f, "test_group_5"), "path_that_does_not_exist", 3.0) == 3.0
@test get!(f["test_group_5"], "path_created_by_get!", 4.0) == 4.0
@test f["test_group_5"]["path_created_by_get!"] == 4.0

# Test that the default element constructor is called only when necessary
# First, we'll access a group which *does* exist, so the default constructor is not called:
@test get(() -> error("this should not be called because the key exists"), f["test_group_5"], "path_created_by_get!") == 4.0

# Then we'll access a few keys that do not exist, and verify that the default
# constructor was called:
default_called = false
get(f, "path_that_does_not_exist") do
    global default_called
    default_called = true
end
@test default_called

default_called = false
get!(JLD2.Group(f, "test_group_6"), "path_created_by_get!") do
    global default_called
    default_called = true
end
@test default_called
close(f)

f = jldopen(fn, "r")
@test read(f, "/test_group_1/x1") == 1
@test read(f, "/test_group_1/x2") == 2
@test read(f, "/test_group_2/x1") == 3
@test read(f, "/test_group_3/x1") == 4
@test read(f, "/test_group_3/contained_group/x1") == 5
@test_throws KeyError read(f, "/test_group_3/contained_group/x2")
@test_throws KeyError read(f, "/test_group_3/contained_group/x2/oops")

@test read(f, "test_group_1/x1") == 1
@test read(f, "test_group_1/x2") == 2
@test read(f, "test_group_2/x1") == 3
@test read(f, "test_group_3/x1") == 4
@test read(f, "test_group_3/contained_group/x1") == 5

@test read(f, "test_group_1")["x1"]  == 1
@test read(f, "test_group_1")["x2"]  == 2
@test read(f, "test_group_2")["x1"]  == 3
@test read(f, "test_group_3")["x1"]  == 4
@test read(f, "test_group_3")["contained_group"]["x1"]  == 5

@test read(f, "test_group_1")["/test_group_2/x1"]  == 3
@test read(f, "test_group_3/contained_group")["x1"]  == 5
@test read(f, "test_group_3")["contained_group/x1"]  == 5

@test !isempty(read(f, "group_with_one_group_child"))
@test !isempty(read(f, "test_group_1"))
@test isempty(read(f, "empty_group"))
@test isempty(read(f, "group_with_one_group_child/empty_group"))

@test get(f, "path_that_does_not_exist", 1.0) == 1.0
@test get(f["test_group_5"], "path_that_does_not_exist", 1.0) == 1.0
@test f["test_group_5"]["path_created_by_get!"] == 4.0

# Throws an ArgumentError because the element doesn't exist and cannot
# be created for a read-only file:
@test_throws ArgumentError get!(f, "path_that_does_not_exist", 1.0)
@test_throws ArgumentError get!(f["test_group_5"], "path_that_does_not_exist", 1.0)
# Does not throw an ArgumentError because the element does exist:
@test get!(f["test_group_1"], "x1", 2) == 1

# Make sure printing doesn't error
show(IOBuffer(), f)
close(f)

# Test for non-default group size
@testset "Custom estimated Group size" begin
    fn = joinpath(mktempdir(), "test.jld2")

    # It is difficult to test the behaviour
    # of non-default group allocation size
    # as results should always be independent of this.
    # Test against regressions by comparing to current number of allocated bytes.
    # This may need to be updated if the relevant code sections are changed.
    # Also test for persistence of the settings when loading a file again.

    # An empty group leaves has 151 bytes by default
    jldopen(fn, "w") do f
        JLD2.Group(f, "g")
    end
    jldopen(fn, "r") do f
        g = f["g"]
        chunk_length = g.last_chunk_checksum_offset - g.last_chunk_start_offset
        @test chunk_length == 151
    end

    # An empty group with non-default size allocates more space
    jldopen(fn, "w") do f
        g = JLD2.Group(f, "g"; est_num_entries=42, est_link_name_len=21)
        @test g.est_link_name_len == 21
        @test g.est_num_entries == 42
    end
    jldopen(fn, "r") do f
        g = f["g"]
        @test g.est_link_name_len == 21
        @test g.est_num_entries == 42
        chunk_length = g.last_chunk_checksum_offset - g.last_chunk_start_offset
        @test chunk_length == 1614
    end
end

# Test for Issue #307
@testset "Repeatedly appending to file" begin
    fn = joinpath(mktempdir(), "test.jld2")
    for n = 1:25
        jldopen(fn, "a") do fid
            group = JLD2.Group(fid, string(n))
            group["result"] = 1
            @test keys(fid) == string.(1:n)
        end
    end
end


text:
struct ForwardDiffFactorization{T,D,N,M,F} <: Factorization{D}
    value::M
    partials::NTuple{N,M}
    factorization::F
end

ForwardDiffFactorization(A) = ForwardDiffFactorization(factorize, A)

function ForwardDiffFactorization(f, A::AbstractMatrix{ForwardDiff.Dual{T,V,N}}) where {T,V,N}
    Av = ForwardDiff.value.(A)
    Ap = ntuple(i->ForwardDiff.partials.(A, i), N)
    Afac = f(Av)
    return ForwardDiffFactorization{T,ForwardDiff.Dual{T,V,N}}(Av, Ap, Afac)
end

function ForwardDiffFactorization{T,D}(Av::M, Ap::NTuple{N,M}, Afac::F) where {T,D,N,M,F}

    return ForwardDiffFactorization{T,D,N,M,F}(Av, Ap, Afac)
end

function Base.adjoint(A::ForwardDiffFactorization{T,D,N,M,F}) where {T,D,N,M,F}
    # value, partials, and factorization of A
    Av = adjoint(A.value)
    Ap = adjoint.(A.partials)
    Afac = adjoint(A.factorization)
    return ForwardDiffFactorization{T,D,N,M,F}(Av, Ap, Afac)
end

function Base.transpose(A::ForwardDiffFactorization{T,D,N,M,F}) where {T,D,N,M,F}
    # value, partials, and factorization of A
    Av = transpose(A.value)
    Ap = transpose.(A.partials)
    Afac = transpose(A.factorization)
    return ForwardDiffFactorization{T,D,N,M,F}(Av, Ap, Afac)
end

function LinearAlgebra.ldiv!(A::ForwardDiffFactorization{T,D,N,M,F}, B) where {T,D,N,M,F}

    # value, partials, and factorization of A
    Av = A.value
    Ap = A.partials
    Afac = A.factorization

    # values and partials of B
    Bv = ForwardDiff.value.(B)
    Bp = ntuple(i->ForwardDiff.partials.(B, i), N)

    # calculate output values, store in Bv
    ldiv!(Afac, Bv)

    # calculate output partials, store in Bp
    ldiv!.(Ref(Afac), mul!.(Bp, Ap, Ref(Bv), -1, 1))

    # overwrite B with result
    for j in eachindex(B)
        B[j] = ForwardDiff.Dual{T}(Bv[j], ntuple(i->Bp[i][j], N)...)
    end

    return B
end

safe_lu(A; check::Bool = true) = lu(A; check)

function safe_lu(A::SparseMatrixCSC{<:ForwardDiff.Dual, <:Any}; check::Bool = true)
    return ForwardDiffFactorization(lu, A)
end


text:
module Shell

export @esc_cmd

SHELL   = Sys.iswindows() ? "cmd" :
          Sys.isapple()   ? get(ENV, "SHELL", "bash") :
          readchomp(pipeline(`getent passwd $(get(ENV, "LOGNAME", "root"))`, `cut -d: -f7`))
CHOMP   = true
SOURCE  = false
CAPTURE = false
DRYRUN  = false
TMPDIR  = tempdir()

"""
    run(cmd::AbstractString; shell=SHELL, capture=CAPTURE, chomp=CHOMP,
        dryrun=DRYRUN, source=SOURCE)

Run your command string in shell.

# Examples
```jldoctest
julia> using Shell

julia> Shell.run(raw"echo \$SHELL", capture=true)
"/bin/zsh"

julia> Shell.run(raw"for i in bu fan; do echo \$i; done")
bu
fan

julia> files = ["temp file 1", "temp file 2"]
2-element Array{String,1}:
 "temp file 1"
 "temp file 2"

julia> filelist = esc`\$files.txt`
"'temp file 1.txt' 'temp file 2.txt'"

julia> Shell.run("touch \$filelist")

julia> Shell.run("ls > 'temp file 0.txt'")

julia> Shell.run("cat 'temp file 0.txt' | grep temp")
temp file 0.txt
temp file 1.txt
temp file 2.txt

julia> Shell.run("rm 'temp file'*")
```

* You should properly escape all special characters manually.
* use `dryrun=true` to check the command to be run without actually running.
* To capture output, set `capture=true`.
* To avoid escaping `\$` everytime, you can use raw string,
  like `raw"echo \$PATH"`
* You can change the default shell (`zsh` in linux and `cmd` in windows)
  using `setshell("other_shell")`.
* In Windows, shell should be `cmd` or `powershell`.
"""
function run(cmd::AbstractString; shell=SHELL, capture=CAPTURE,
             chomp=CHOMP, dryrun=DRYRUN, source=SOURCE)
    result = nothing
    file = joinpath(TMPDIR, basename(tempname()))
    pre_script = ""
    command = ``
    if Sys.iswindows()
        if shell == "cmd"
            file = "$file.bat"
            pre_script = "@echo off\nchcp 65001 >nul"
            command = `$file`
        elseif shell == "powershell"
            file = "$file.ps1"
            pre_script = raw"chcp 65001 >$null"
            command = `powershell -command $file`
        elseif shell == "wsl"
            wslshell = split(readchomp(`wsl getent passwd \$LOGNAME`), ":")[end]
            file = "$file.sh"
            source && (pre_script = shell_source(wslshell))
            ffile = replace(file, "\\" => "\\\\")
            wslfile = readchomp(`wsl wslpath $ffile`)
            command = `wsl $wslshell $wslfile`
        else
            @error("Only support cmd/powershell/wsl in Windows.")
        end
    else
        source && (pre_script = shell_source(shell))
        command = `$shell $file`
    end

    open(file, "w") do f
        println(f, pre_script)
        println(f, cmd)
        flush(f)
    end

    if dryrun
        result = println(strip(read(file, String)))
    elseif capture
        result = chomp ? readchomp(command) : read(command, String)
    else
        Base.run(command)
    end
    rm(file)
    return result
end

function runfile(file; background=true, shell=SHELL)
    if background
        return Base.run(detach(`$shell $file`), wait=false)
    else
        return Base.run(`$shell $file`)
    end
end

shell_source(shell) =
    (endswith(shell, "/zsh") || shell == "zsh") ? "[ -f ~/.zshrc ] && source ~/.zshrc >/dev/null" :
    (endswith(shell, "/bash") || shell == "bash") ? "[ -f ~/.bashrc ] && source ~/.bashrc" :
    (endswith(shell, "/fish") || shell == "fish") ? "[ -f ~/.config/fish/config.fish ] && source ~/.config/fish/config.fish" :
    (@warn("Please make a PR to support source your awesome `$shell` shell!"); "")

"""
    setshell(shell::AbstractString)

Specify which shell to use (Windows defaults to `cmd`, and `zsh` otherwise).
"""
function setshell(shell::AbstractString)
    global SHELL
    SHELL = shell
end

"""
    settmp(tmp::AbstractString)

Set the temp folder, defaults to `/tmp/`.
"""
function settmp(tmp::AbstractString)
    global TMPDIR
    TMPDIR = tmp
    mkpath(TMPDIR)
end

"""
    setchomp(chomp::Bool)

Set whether chomp the output (default is true).
"""
function setchomp(chomp::Bool)
    global CHOMP
    CHOMP = chomp
end

"""
    setissource(source::Bool)

Whether source the rc file (e.g. `.zshrc`) before run script (default is true).
"""
function setissource(source::Bool)
    global SOURCE
    SOURCE = source
end

"""
    setisdryrun(dryrun::Bool)

Whether dryrun the rc file (e.g. `.zshrc`) before run script (default is true).
"""
function setisdryrun(dryrun::Bool)
    global DRYRUN
    DRYRUN = dryrun
end

"""
    setiscapture(capture::Bool)

Whether capture the rc file (e.g. `.zshrc`) before run script (default is true).
"""
function setiscapture(capture::Bool)
    global CAPTURE
    CAPTURE = capture
end

"""
    @esc_str -> String

Help you escape special characters for the shell.

# Examples
```jldoctest
julia> files = ["temp file 1", "temp file 2"]
2-element Array{String,1}:
 "temp file 1"
 "temp file 2"

julia> filelist = esc`\$files.txt`
"'temp file 1.txt' 'temp file 2.txt'"

julia> Shell.run("touch \$filelist")

julia> Shell.run("rm \$filelist")
```

Be careful, the escape treat space separated terms individually.
Put them into a varible to get properly escaped.

# Examples
```jldoctest
julia> esc`temp file 0.txt`
"temp file 0.txt"

julia> file = "temp file 0.txt"
"temp file 0.txt"

julia> esc`\$file`
"'temp file 0.txt'"
```

* Not working for `cmd` in Windows because it treats single quotes differently.
"""
macro esc_cmd(cmd)
    esc(:(join(map((@cmd $cmd).exec) do arg
        replace(sprint() do io
            Base.print_shell_word(io, arg, Base.shell_special)
        end, '`' => "\\`")
    end, ' ')))
end

end # module


text:
"""
resample(uvals::Vector{AbstractUncertainValue})

Treat `uvals` as a dataset and resample each value of `uvals` once.
Returns an `length(uvals)`-element vector.
"""
resample(uvals::Vector{AbstractUncertainValue}) = resample.(uvals)

"""
resample(uvals::Vector{AbstractUncertainValue}, n::Int)

Treat `uvals` as a dataset and resample it `n` times.
Returns `n` resampled draws of `uvals`, each being a `length(uvals)`-element vector.
For each returned vector, the i-th element is a unique draw of `uvals[i]`.
"""
function resample(uvals::Vector{AbstractUncertainValue}, n::Int)
    [resample.(uvals) for i = 1:n]
end

"""
resample_elwise(uvals::Vector{AbstractUncertainValue}, n::Int)

Resample each element in `uvals` `n` times. The i-th entry in the returned
vector is a `n`-element vector consisting of `n` unique draws of `uvals[i]`.
"""
function resample_elwise(uvals::Vector{AbstractUncertainValue}, n::Int)
    [resample(uvals[i], n) for i = 1:length(uvals)]
end


text:
using LinearAlgebra

mutable struct EpicomposeQuadratic{F, M, P, V, R} <: Epicompose
    L::M
    Q::P
    q::V
    gamma::Maybe{R}
    fact::Maybe{F}
    function EpicomposeQuadratic{F, M, P, V, R}(L::M, Q::P, q::V) where {
        R <: Real,
        M <: AbstractMatrix{R},
        P <: AbstractMatrix{R},
        V <: AbstractVector{R},
        F <: Factorization,
    }
        new(L, Q, q, nothing, nothing)
    end
end

function EpicomposeQuadratic(L, Q::P, q) where {R, P <: DenseMatrix{R}}
    return EpicomposeQuadratic{
        LinearAlgebra.Cholesky{R}, typeof(L), P, typeof(q), real(eltype(L))
    }(L, Q, q)
end

function EpicomposeQuadratic(L, Q::P, q) where {R, P <: SparseMatrixCSC{R}}
    return EpicomposeQuadratic{
        SuiteSparse.CHOLMOD.Factor{R}, typeof(L), P, typeof(q), real(eltype(L))
    }(L, Q, q)
end

# TODO: enable construction from other types of quadratics, e.g. LeastSquares
# TODO: probably some access methods are needed to obtain Hessian and linear term?

function EpicomposeQuadratic(L, f::Q) where {Q <: Quadratic}
    return EpicomposeQuadratic(L, f.Q, f.q)
end

function factor_step!(g::EpicomposeQuadratic, gamma::R) where {R <: Real}
    g.gamma = gamma
    g.fact = cholesky(g.Q + (g.L' * g.L)/gamma)
end

function prox!(y::AbstractArray{R}, g::EpicomposeQuadratic, x::AbstractArray{R}, gamma::R=one(R)) where {
    R <: Real
}
    if g.gamma === nothing || !isapprox(gamma, g.gamma)
        factor_step!(g, gamma)
    end
    p = g.fact\((g.L' * x) / gamma - g.q)
    fy = dot(p, g.Q * p)/2 + dot(p, g.q)
    mul!(y, g.L, p)
    return fy
end

function prox_naive(g::EpicomposeQuadratic, x::AbstractArray{R}, gamma::R=one(R)) where {
    R <: Real
}
    S = g.Q + (g.L' * g.L) / gamma
    p = S\((g.L' * x) / gamma - g.q)
    fy = dot(p, g.Q * p)/2 + dot(p, g.q)
    y = g.L * p
    return y, fy
end


text:
# Welcome.

# This is a regular file, which you can edit as much as you want,
# then save for later or discard. You can open a fresh copy by
# pressing `Ctrl-Space` and typing `tutorial` followed by `Enter`.

# Use the "toggle console" command to see any output. If you have
# issues, please let me know at
# http://discuss.junolab.org/

# ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––

# In Juno we'll be using the Julia language. This is a Julia file and
# lines beginning with `#` are comments.

# Here's a very simple code block:

3 * 6 + 4

# We can evaluate code by right clicking on it and selecting "Evaluate".
# You should see the result, "22", pop up next to the code.

# It's quicker to use Ctrl-Enter (or Cmd-Enter on Mac). Try adding "+2"
# to this line and then evaluating.

2+2

# We can also define our own functions, such as this one which doubles a
# number:

function double(x)
  return 2x
end

# Click anywhere on the function definition and press `Ctrl-Enter`
# again. You should see a tick pop up next to the "end" to show that
# the function has been defined.

# Now we can call our function:

double(-4)

# After evaluating this you should see the result "20" pop up. You may
# also notice that the "10" becomes underlined. Try clicking and
# dragging it from side to side to double another number.

# We can do more advanced things, like dealing with matrices. Try

x = rand(5, 4)

# If you evaluate this you'll see "Matrix Float64, 5×5" pop up next to it.
# A bold header like this indicates that you can click on the box
# to expand it and see the whole matrix.

# If you're not sure what a function does it's easy to find out. Just right
# click on the function's name (like "rand" above) and click "Toggle Docs".
# The shortcut for this is Ctrl-D. You can also remove the documentation
# by right clicking and selecting "Remove result".

# Let's try installing a package, then do some visualisation.

Pkg.add("Images")

# If you start to type `Images` between the quotes, you should see the
# autocomplete pop up with a list of packages. Select `Images` with tab
# and evaluate to install the package. You'll see a working indicator
# in the bottom-left corner as the package is loaded. Type `Ctrl-Space`
# followed `toggle console` and `Enter` to see the output.

# Now you can load the package.

using Images

# Let's try visualising this function. You can type latex symbols by
# typing a `\` and selecting from the autocomplete, e.g. `\phi` or `\le`.
τ = 2π
# You can use `Shift-Enter` to quickly evaluate definitions.

const ϕ = golden

function foo(z)
  c = (φ-2)+(φ-1)im
  max = 80
  for n = 1:max
    abs(z) ≥ 2 && return n-1
    z = z^2 + c
  end
  return max
end

# Let's try this out – play around with the numbers!

foo(-1)
foo(5+0.75im)

# Let's apply foo to a grid of numbers. Don't worry too much if you don't
# understand this code, but make sure to evaluate it. (Bonus points for
# using that Ctrl-D trick here, though)

foo(x, y) = foo(x + y*im)

foo_grid(n) =
  broadcast(foo,
            linspace(-0.5, 1, n)',
            linspace(-1, 0.5, n))

# Can you see the pattern?

foo_grid(10)

# Let's try it as an image:

convert(Image, scale(foo_grid(500), 1/80))

# It's a Julia set, of course!

# If you expand the image you can also click on the `Dictionary`
# field to see some interesting properties about it (if that's the
# kind of thing that floats your goat).

# If you want to learn more about Julia, take a look at
# http://julialang.org/learning/

# There is also the Julia manual and the Juno documentation, both
# of which are available from the help menu above.

# Enjoy!