text:
function plot3(x,y,z,lineSpec::String;options=Dict())
    MATLAB.@mput x
    MATLAB.@mput y
    MATLAB.@mput z
    optionsString = namevaluepairs(options)
    plotString = "plot3(x,y,z"*",\"$lineSpec\""*optionsString*")"
    MATLAB.eval_string(plotString)
end

function plot3(x,y,z;options=Dict())
    MATLAB.@mput x
    MATLAB.@mput y
    MATLAB.@mput z
    optionsString = namevaluepairs(options)
    plotString = "plot3(x,y,z"*optionsString*")"
    MATLAB.eval_string(plotString)
end

#TODO: Improve 3D Plotting
function surf(X,Y,Z;options=Dict())
    MATLAB.@mput X
    MATLAB.@mput Y
    MATLAB.@mput Z
    optionsString = namevaluepairs(options)
    plotString = "surf(X,Y,Z"*optionsString*")"
    MATLAB.eval_string(plotString)
end

function surf(X,Y,Z,lineSpec;options=Dict())
    MATLAB.@mput X
    MATLAB.@mput Y
    MATLAB.@mput Z
    optionsString = namevaluepairs(options)
    plotString = "surf(X,Y,Z"*",\"$lineSpec\""*optionsString*")"
    MATLAB.eval_string(plotString)
end

function surface(X,Y,Z;options=Dict())
    MATLAB.@mput X
    MATLAB.@mput Y
    MATLAB.@mput Z
    optionsString = namevaluepairs(options)
    plotString = "surface(X,Y,Z"*optionsString*")"
    MATLAB.eval_string(plotString)
end

function surface(X,Y,Z,lineSpec;options=Dict())
    MATLAB.@mput X
    MATLAB.@mput Y
    MATLAB.@mput Z
    optionsString = namevaluepairs(options)
    plotString = "surface(X,Y,Z"*",\"$lineSpec\""*optionsString*")"
    MATLAB.eval_string(plotString)
end

function contourf(X,Y,Z;options=Dict())
    MATLAB.@mput X
    MATLAB.@mput Y
    MATLAB.@mput Z
    optionsString = namevaluepairs(options)
    plotString = "contourf(X,Y,Z"*optionsString*")"
    MATLAB.eval_string(plotString)
end

function contourf(X,Y,Z,lineSpec;options=Dict())
    MATLAB.@mput X
    MATLAB.@mput Y
    MATLAB.@mput Z
    optionsString = namevaluepairs(options)
    plotString = "contourf(X,Y,Z"*",\"$lineSpec\""*optionsString*")"
    MATLAB.eval_string(plotString)
end

function contour(X,Y,Z;options=Dict())
    MATLAB.@mput X
    MATLAB.@mput Y
    MATLAB.@mput Z
    optionsString = namevaluepairs(options)
    plotString = "contour(X,Y,Z"*optionsString*")"
    MATLAB.eval_string(plotString)
end

function contour(X,Y,Z,lineSpec;options=Dict())
    MATLAB.@mput X
    MATLAB.@mput Y
    MATLAB.@mput Z
    optionsString = namevaluepairs(options)
    plotString = "contour(X,Y,Z"*",\"$lineSpec\""*optionsString*")"
    MATLAB.eval_string(plotString)
end


text:
module SimpleHeatmaps

export
    current_colormap,
    set_colormap,
    matshow,
    openanim,
    addframe

include("matshow.jl")

end # module

ext: jl


text:
using Distributions

module ConjugatePriors
	export update
	using Distributions

	include("bayes.jl")
end


text:
s = read_SSR(   joinpath(dirname(@__FILE__), "..", "data", "test_Hz19.5-testing.bdf"))
s = read_evt(s, joinpath(dirname(@__FILE__), "..", "data", "test.evt"))


text:
mutable struct Drive{V} <: State{V}
    value::V
    array::Vector{V}
    tick::Int
end

Drive(; tick, unit, _value, _type, _...) = begin
    t = value(tick)
    U = value(unit)
    V = valuetype(_type, U)
    a = _value
    v = a[1]
    Drive{V}(v, a, t)
end

constructortags(::Val{:Drive}) = (:tick, :unit)

updatetags!(d, ::Val{:Drive}; _...) = begin
    !haskey(d, :tick) && (d[:tick] = :(context.clock.tick))
end

genvartype(v::VarInfo, ::Val{:Drive}; V, _...) = @q Drive{$V}

gendefault(v::VarInfo, ::Val{:Drive}) = begin
    s = gettag(v, :from)
    x = if isnothing(s)
        k = gettag(v, :by)
        !isnothing(k) && error("missing `from` provider for `by` = $k")
        istag(v, :parameter) ? genparameter(v) : genbody(v)
    else
        istag(v, :parameter) && error("`parameter` is not allowed with provider: $s")
        !isnothing(v.body) && error("function body is not allowed with provider: $s\n$(v.body)")
        #HACK: needs quot if key is a symbol from VarInfo name
        k = gettag(v, :by, Meta.quot(v.name))
        @q $C.value($s)[!, $k]
    end
    u = gettag(v, :unit)
    @q $C.unitfy($x, $C.value($u))
end

genupdate(v::VarInfo, ::Val{:Drive}, ::MainStep; kw...) = begin
    t = gettag(v, :tick)
    @gensym s t0 t1 i e
    @q let $s = $(symstate(v)),
           $t0 = $s.tick,
           $t1 = $C.value($t),
           $i = $t1 - $t0 + 1,
           $e = $s.array[$i]
        $C.store!($s, $e)
    end
end


text:
################################
###    Summary statistics    ###
################################
"""
	poisson_fixation(empirical_values,λds, λdn)

Divergence sampling from Poisson distribution. The expected neutral and selected fixations are subset through their relative expected rates ([`fixNeut`](@ref), [`fixNegB`](@ref), [`fixPosSim`](@ref)). Empirical values are used are used to simulate the locus *L* along a branch of time *T* from which the expected *Ds* and *Dn* raw count estimated given the mutation rate (``\\mu``). Random number generation is used to subset samples arbitrarily given the success rate ``\\lambda`` in the distribution.

```math
\\mathbb{E}[D_N] = X \\in Poisson\\left(\\lambda = D \\times \\left[\\frac{\\mathbb{E}[D_+] + \\mathbb{E}[D_-]}{\\mathbb{E}[D_+] + \\mathbb{E}[D_-] + \\mathbb{E}[D_0]}\\right]\\right)
```
```math
\\mathbb{E}[D_S] = X \\in Poisson\\left(\\lambda = D \\times \\left[\\frac{\\mathbb{E}[D_0]}{\\mathbb{E}[D_+] + \\mathbb{E}[D_-] + \\mathbb{E}[D_0]}\\right]\\right)
```
# Arguments
 - `observedValues::Array`: Array containing the total observed divergence.
 - ` λds::Float64`: expected neutral fixations rate.
 - ` λdn::Float64`: expected selected fixations rate.
# Returns
 - `Array{Int64,1}` containing the expected count of neutral and selected fixations.

"""
function poisson_fixation(;empirical_values::Array, λds::Array, λdn::Array,λweak::Array,λstrong::Array)

    ds             = @. λds / (λds + λdn) * empirical_values
    dn             = @. λdn / (λds + λdn) * empirical_values
    dn_weak        = @. λweak / (λds + λdn) * empirical_values
    dn_strong      = @. λstrong / (λds + λdn) * empirical_values

    sampled_Ds     = pois_rand.(ds)
    sampled_Dn     = pois_rand.(dn)
    sampled_weak   = pois_rand.(dn_weak)
    sampled_strong = pois_rand.(dn_strong)

    alphas         = @. [sampled_weak/sampled_Dn sampled_strong/sampled_Dn (sampled_weak+sampled_strong)/sampled_Dn]

	out            = alphas,sampled_Dn, sampled_Ds
	return out
end

"""
	poisson_polymorphism(empirical_values,λps,λpn)

Polymorphism sampling from Poisson distributions. The total expected neutral and selected polimorphism are subset through the relative expected rates at the frequency spectrum ([`fixNeut`](@ref), [`DiscSFSNeutDown`](@ref),). Empirical sfs are used to simulate the locus *L* along a branch of time *T* from which the expected *Ps* and *Pn* raw count are estimated given the mutation rate (``\\mu``). Random number generation is used to subset samples arbitrarily from the whole sfs given each frequency success rate ``\\lambda`` in the distribution.

The success rate managing the Poisson distribution by the observed count each frequency.  We considered both sampling variance and process variance is affecting the number of variable alleles we sample from SFS. This variance arises from the random mutation-fixation process along the branch. To incorporate this variance we do one sample per frequency-bin and use the total sampled variation and the SFS for the summary statistics.

```math
\\mathbb{E}[P_N] = \\sum_{x=0}^{x=1} X \\in Poisson\\left(\\lambda = SFS_{(x)} \\times \\left[\\frac{\\mathbb{E}[P_{+(x)}] + \\mathbb{E}[P_{-(x)}]}{\\mathbb{E}[P_{+(x)}] + \\mathbb{E}[P_{-(x)}] + \\mathbb{E}[P_{0(x)}]}\\right]\\right)
```

```math
\\mathbb{E}[P_S] = \\sum_{x=0}^{x=1} X \\in Poisson\\left(\\lambda = SFS_{(x)} \\times \\left[\\frac{\\mathbb{E}[P_{0(x)}]}{\\mathbb{E}[P_{+(x)}] + \\mathbb{E}[P_{-(x)}] + \\mathbb{E}[P_{0(x)}]}\\right]\\right)
```

# Arguments
 - `observedValues::Array{Int64,1}`: Array containing the total observed divergence.
 - ` λps::Array{Float64,1} `: expected neutral site frequency spectrum rate.
 - ` λpn::Array{Float64,1} `: expected selected site frequency spectrum rate.
# Returns
 - `Array{Int64,2}` containing the expected total count of neutral and selected polymorphism.

"""
function poisson_polymorphism(;empirical_values::Array, λps::Array, λpn::Array,pol::Array)

	# Neutral λ;156948
	λ1 = @. λps / (λps + λpn) * empirical_values'
	# Selected λ;
	λ2 = @. λpn / (λps + λpn) * empirical_values'

	# Neutral λ;
	#=λ1 = @. λps / pol[:,1] * observedValues[:,2]'
	# Selected λ;
	λ2 = @. λpn / pol[:,2] * observedValues[:,1]'=#

	# Replace negative and 0 rates values. Strange behaviour depending on θ and Γ values
	# Relative rates output NaN values due to 0 divisons.

	replace!(λ1,NaN=>1)
	replace!(λ2,NaN=>1)

	sampled_Ps = pois_rand.(λ1)
	sampled_Pn = pois_rand.(λ2)

	return (sampled_Pn, sampled_Ps)
end

"""

	Ouput the expected values from the Poisson sampling process. Please check [`poisson_fixation`](@ref) and [`poisson_polymorphism`](@ref) to understand the samplingn process. α(x) is estimated through the expected values of Dn, Ds, Pn and Ps.

	# Arguments
	 - `afs::Vector`: Array containing SFS subset
	 - `divergence::Vector`: Array containing divergence counts
	 - `observed::Matrix`: Array containing SFS
	 - `m::Array`: Array containing modeled parameters
	 - `nt::Array`: Array containing neutral polymorphic rates at frequencies dac
	 - `sl::Array`: Array containing selective polymorphic rates at frequencies dac
	 - `d::Array`: Array containing fixation rates
	 - `p::Array`: Array containing total polymorphic rates
	 - `dac::Array`: Array containing Derived Allel Counts to subset
	# Returns
	Expected α_x values
	 - `Array{Int64,2}` true α, γ, shape parameter and α_x values.

	sampling_from_rates(afs,divergence,observed,m,nt,d,p,dac)
"""
function sampling_from_rates(afs::Vector,divergence::Vector,observed::Matrix,m::Array,nt::Array,sl::Array,d::Array,p::Array,dac::Array)

	ds             = d[:,1]
	dn             = d[:,2]
	dn_weak        = d[:,3]
	dn_strong      = d[:,4]
	gn             = abs.(m[:,4])
	sh             = round.(m[:,end-1],digits=5)

	#alxSummStat, alphasDiv, expectedDn, expectedDs, expectedPn, expectedPs = sampledAlpha(divergence=divergence,afs=s,observed=observed,λdiv=[ds,dn,dn_weak,dn_strong],λpol=[permutedims(nt),permutedims(sl)],dac=dnac)

	## Outputs
	alphas, expDn, expDs = poisson_fixation(empirical_values=divergence,λds=ds,λdn=dn,λweak=dn_weak,λstrong=dn_strong)
	expPn, expPs         = poisson_polymorphism(empirical_values=afs,λps=nt,λpn=sl,pol=p)

	θ_pn = @. observed[1,1] / (expPn[:,1]);
	θ_ps = @. observed[1,2] / (expPs[:,1]);

	r_pn = ones(size(expPn));
	r_ps = ones(size(expPs));
	for j=2:size(dac,1)
		@.  r_pn[:,j] = (observed[j,1]) / (expPn[:,j]);
		@. r_pn[:,j] = r_pn[:,j] / θ_pn;

		@. r_ps[:,j] = (observed[j,2]) / (expPs[:,j]);
		@. r_ps[:,j] = r_ps[:,j] / θ_ps;
	end

	expPn_rj = zeros(size(expPn))
	expPs_rj = zeros(size(expPs))

	expPn_rj[:,1] = @. expPn[:,1] * θ_pn;
	expPs_rj[:,1] = @. expPs[:,1] * θ_ps;
	for j = 2:size(dac,1)
		@. expPn_rj[:,j] =  expPn[:,j] * θ_pn * r_pn[:,j];
		@. expPs_rj[:,j] =  expPs[:,j] * θ_ps * r_ps[:,j];
	end

	## Alpha from expected values. Used as summary statistics
    α_x    = @. round(1 - ((expDs/expDn) * (expPn/expPs)),digits=5)
    α_x_rj	 = @. round(1 - ((expDs/expDn) * (expPn_rj/expPs_rj)),digits=5)

	expected_values = hcat(round.(alphas,digits=5),gn,sh,α_x)

	return(expected_values)
end

"""
	summary_statistics(param::parameters,rates::JLD2.JLDFile,analysis_folder::String,summstat_size::Int64,replicas::Int64,bootstrap::Bool)

Estimate summary statistics using observed data and analytical rates. *analysis_folder* will check for the SFS and divergence file and will be used to output summary statistics

# Arguments
 - `param::parameters` : Mutable structure containing the models
 - `rates::JLD2.JLDFile` : HDF5 containing solved models, fixation and polymorphic rates
 - `analysis_folder::String` : Folder containing the SFS and divergence files. It will be used to output the observed data and summary estatistics.
 - `summstat_size::Int64` : Number of summary statistics
 - `replicas::Int64` : Number of bootstrap replicas
 - `bootstrap::Bool` : Boolean to perform or not bootstrapping
# Output
 - Obserded data and summary statistics to ABC inference
"""
function summary_statistics(;param::parameters,h5_file::String,analysis_folder::String,summstat_size::Int64,replicas::Int64=1,bootstrap::Bool=false)

	# param = parameters(n=50,dac=[1,2,4,5,10,20,25,30,50,70]);
	# analysis_folder="/home/jmurga/test/simulations/prueba/";
	# # h5_file="/home/jmurga/test/simulations/prueba/prueba.jld2";
	# h5_file="/home/jmurga/test/simulations/rates.jld2";
	# replicas=1;
	# bootstrap=true;
	# summstat_size=10^5

	#Opening files
	sfs_files        = filter(x -> occursin("sfs",x), readdir(analysis_folder,join=true));
	divergence_files = filter(x -> occursin("div",x), readdir(analysis_folder,join=true));

	sfs,divergence,α,observed = open_sfs_div(sfs_files,divergence_files,param.dac,replicas,bootstrap);

	#Open rates
	h5    = jldopen(h5_file);
	h5_subset   = h5[string(param.N) * "/" * string(param.n)];

	#Subset index
	idx    = sample(1:size(h5_subset["models"],1),summstat_size,replace=false);
	models = Array(view(h5_subset["models"],idx,:));
	dsdn   = Array(view(h5_subset["dsdn"],idx,:));
	pol    = Array(view(h5_subset["pol"],idx,:));

	# Filtering polymorphic rate by dac
	n    = hcat(map(x -> view(h5_subset["neut"][x],:),param.dac)...);
	s    = hcat(map(x -> view(h5_subset["sel"][x],:),param.dac)...);
	neut = Array(view(n,idx,:));
	sel  = Array(view(s,idx,:));

	#Making summaries
	f(afs,divergence,observed,m=models,nt=neut,sl=sel,d=dsdn,p=pol,dac=param.dac) = sampling_from_rates(afs,divergence,observed,m,nt,sl,d,p,dac);

	expected_values = [hcat(models[:,1:5],@. round(1 - (dsdn[:,1]/dsdn[:,2] * sel/neut),digits=5))]

	expected_values = progress_pmap(f,sfs,divergence,observed;progress = Progress(size(sfs,1),desc="Sampling αₓ "));

	# Controling outlier cases
	expected_values = pmapbatch(filter_expected,expected_values);

	# nuisance =  [[1.16  0.62  0.32  0.04  1.  1.  1. 1.  1.  1.],[1.08  0.7   0.2   0.41  1.51  1.3   1.18  1.01  0.94  1.04],[1.16  0.62  0.32  0.04  1.11  1.12  1.03  0.92  0.93  1.02],[0.76  0.91  1.44  1.98  0.22  0.46  0.58  0.7  0.77  1.01],[0.8   1.03  1.91  3.02  1.6   0.45  0.16  0.05  0.14  0.26]];

	# expected_values = progress_pmap(y -> vcat(y,vcat(map(x-> nuisance_alpha(y,x),nuisance)...)),expected_values)

	# flt_b(x) = any(x .> 1)
	# flt = pmapbatch(y -> .!mapslices(flt_b,y[:,6:end],dims=2)[:,1],expected_values)

	# flt_v = (x,y) -> @view x[y,:]

	# expected_values = pmapbatch(flt_v,expected_values,flt)

	w(x,name) = CSV.write(name,DataFrame(x,:auto),delim='\t',header=false);

	progress_map(w, α, analysis_folder * "/alphas_" .* string.(1:size(sfs,1)) .* ".txt";progress= Progress(size(sfs_files,1),desc="Writting α "));

	progress_pmap(w, expected_values,  analysis_folder * "/summstat_" .* string.(1:size(sfs,1)) .* ".txt";progress= Progress(size(sfs,1),desc="Writting summaries "));

	return(expected_values)
end

function filter_expected(x::Matrix{Float64})

	replace!(x, -Inf=>NaN)
	x = x[vec(.!any(isnan.(x),dims=2)),:]
	x = x[(x[:,3] .< 1),:]

	return(x)
end

# function nuisance_alpha(a::Matrix{Float64},b::Matrix{Float64})

# 		tmp_param = @view a[:,1:5];
# 		tmp_summs = @view a[:,6:end];

# 		tmp = round.(tmp_summs .* b,digits=5)

# 		out = hcat(tmp_param,tmp);
# 		return(out)
# end


text:
using Base.Dates
using BusinessDays
using DataFrames

bd = BusinessDays

hc_vec = [:BRSettlement, :USSettlement, :USNYSE, :USGovernmentBond, :UKSettlement, :CanadaSettlement, :CanadaTSX, :TARGET]

d0 = Date(1960,01,04)
d1 = Date(2100,01,04)
bd.initcache(hc_vec, d0, d1)

d1vec = collect(d0:d1)
d0vec = fill(d0, length(d1vec))

!isdir("csv") && mkdir("csv")

for hc in hc_vec
	b = convert(Array{Int, 1}, isbday(hc, d1vec))
	writetable(string("csv/julia-isbday-", string(hc), ".csv"), DataFrame(D = d1vec, ISBDAY = b), quotemark = ' ')
end


text:
@testset "Sorts" begin
    x=[3,5,1,4,2]
    BubbleSort!(x)
    @test x == [1,2,3,4,5]
    x=[3,5,1,4,2]
    InsertionSort!(x)
    @test x == [1,2,3,4,5]
    x=[3,5,1,4,2]
    SelectionSort!(x)
    @test x == [1,2,3,4,5]
end

max_stars_repo_name: ashwani-rathee/Julia

max_issues_repo_name: ashwani-rathee/Julia

max_forks_repo_name: ashwani-rathee/Julia


text:
function checkcollision(b::Board)
    x,y = b.location
    dy, dx = size(rotatedtile(b)) .-1
    x+dx <= 10 && x >= 1 && y+dy <= 20 && !any((b[] .!= 0 ) .& (rotatedtile(b) .!= 0))
end

function affine!(b::Board, rotation=0, translation=[0,0])
    oldboard = copy(b)

    b[] -= rotatedtile(b)
    b.orientation -= rotation
    b.location += translation
    if checkcollision(b)
        b[] += rotatedtile(b)
        update_board!(oldboard, b)
        return true
    end
    b.orientation += rotation
    b.location -= translation
    b[] += rotatedtile(b)
    return false
end

rot_right!(b::Board)  = affine!(b, -1, [ 0, 0])
rot_left!(b::Board)   = affine!(b,  1, [ 0, 0])
move_right!(b::Board) = affine!(b,  0, [ 1, 0])
move_left!(b::Board)  = affine!(b,  0, [-1, 0])
drop!(b::Board)       = affine!(b,  0, [ 0, 1])
fast_drop!(b::Board)  = (while drop!(b) end; return false)

function hold!(b::Board)
    if b.allowhold
        oldboard = copy(b)

        b[] -= rotatedtile(b)
        b.tile, b.holdtile = b.holdtile, b.tile
        b.orientation = 0
        b.location = start_location(b.tile)
        b[] += rotatedtile(b)
        print_hold(b)
        b.allowhold = false

        update_board!(oldboard, b)
    end
end

max_stars_repo_name: ChristianKurz/REPLTetris

max_issues_repo_name: ChristianKurz/REPLTetris

max_forks_repo_name: ChristianKurz/REPLTetris


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayProbePropertiesFormat <: SwaggerModel
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol
    host::Any # spec type: Union{ Nothing, String } # spec name: host
    path::Any # spec type: Union{ Nothing, String } # spec name: path
    interval::Any # spec type: Union{ Nothing, Int32 } # spec name: interval
    timeout::Any # spec type: Union{ Nothing, Int32 } # spec name: timeout
    unhealthyThreshold::Any # spec type: Union{ Nothing, Int32 } # spec name: unhealthyThreshold
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function ApplicationGatewayProbePropertiesFormat(;protocol=nothing, host=nothing, path=nothing, interval=nothing, timeout=nothing, unhealthyThreshold=nothing, provisioningState=nothing)
        o = new()
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("host"), host)
        setfield!(o, Symbol("host"), host)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("path"), path)
        setfield!(o, Symbol("path"), path)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("interval"), interval)
        setfield!(o, Symbol("interval"), interval)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("timeout"), timeout)
        setfield!(o, Symbol("timeout"), timeout)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("unhealthyThreshold"), unhealthyThreshold)
        setfield!(o, Symbol("unhealthyThreshold"), unhealthyThreshold)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ApplicationGatewayProbePropertiesFormat

const _property_map_ApplicationGatewayProbePropertiesFormat = Dict{Symbol,Symbol}(Symbol("protocol")=>Symbol("protocol"), Symbol("host")=>Symbol("host"), Symbol("path")=>Symbol("path"), Symbol("interval")=>Symbol("interval"), Symbol("timeout")=>Symbol("timeout"), Symbol("unhealthyThreshold")=>Symbol("unhealthyThreshold"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ApplicationGatewayProbePropertiesFormat = Dict{Symbol,String}(Symbol("protocol")=>"String", Symbol("host")=>"String", Symbol("path")=>"String", Symbol("interval")=>"Int32", Symbol("timeout")=>"Int32", Symbol("unhealthyThreshold")=>"Int32", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ ApplicationGatewayProbePropertiesFormat }) = collect(keys(_property_map_ApplicationGatewayProbePropertiesFormat))
Swagger.property_type(::Type{ ApplicationGatewayProbePropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ApplicationGatewayProbePropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationGatewayProbePropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationGatewayProbePropertiesFormat[property_name]

const _allowed_ApplicationGatewayProbePropertiesFormat_protocol = ["Http", "Https"]

function check_required(o::ApplicationGatewayProbePropertiesFormat)
    true
end

function validate_property(::Type{ ApplicationGatewayProbePropertiesFormat }, name::Symbol, val)
    if name === Symbol("protocol")
        Swagger.validate_param(name, "ApplicationGatewayProbePropertiesFormat", :enum, val, _allowed_ApplicationGatewayProbePropertiesFormat_protocol)
    end
end


text:
# Part of submodule Utils of BetaML _ the Beta Machine Learning Toolkit
# Miscelaneaous funcitons / types


text:
struct ImmutableListNode{T}
    value::T
    tail::Union{ImmutableListNode{T},Nothing}
end

const ImmutableList{T} = Union{ImmutableListNode{T},Nothing}

@inline pushfirst(list, @nospecialize(value)) = ImmutableListNode(value, list)

Base.convert(::Type{ImmutableListNode{T}}, list::ImmutableListNode) where {T} =
    ImmutableListNode{T}(list.value, list.tail)

ilist(xs...) = foldr(ImmutableListNode, xs; init = nothing)

Base.eltype(::Type{ImmutableListNode{T}}) where {T} = T
Base.IteratorSize(::Type{ImmutableListNode{T}}) where {T} = Base.SizeUnknown()
Base.iterate(list::ImmutableListNode{T}, node = list) where {T} =
    if node === nothing
        node
    else
        (node.value, node.tail)
    end

pushsortedby(_, nothing, value) = ImmutableListNode(value, nothing)
function pushsortedby(f, list::ImmutableListNode{T}, value) where {T}
    if isless(f(value), f(list.value))
        return ImmutableListNode{T}(value, list)
    else
        return ImmutableListNode{T}(list.value, pushsortedby(f, list.tail, value))
    end
end

combinesortedby(_, ::Nothing, ::Nothing) = nothing
combinesortedby(_, xs::ImmutableListNode, ::Nothing) = xs
combinesortedby(_, ::Nothing, ys::ImmutableListNode) = ys
function combinesortedby(f, xs::ImmutableListNode{T}, ys::ImmutableListNode{T}) where {T}
    if isless(f(xs.value), f(ys.value))
        return ImmutableListNode{T}(xs.value, combinesortedby(f, xs.tail, ys))
    else
        return ImmutableListNode{T}(ys.value, combinesortedby(f, xs, ys.tail))
    end
end

combine(::Nothing, ::Nothing) = nothing
combine(list::ImmutableListNode, ::Nothing) = list
combine(::Nothing, list::ImmutableListNode) = list
combine(xs::ImmutableListNode{T}, ys::ImmutableListNode{T}) where {T} =
    ImmutableListNode{T}(xs.value, combine(xs.tail, ys))

has(_list, _x) = false
has(list::ImmutableListNode{T}, x::T) where {T} = list.value === x || has(list.tail, x)


text:
using BinDeps
include("./cpuarch_detect.jl")

@BinDeps.setup

sleef = library_dependency("sleef")

provides(Sources, URI("http://shibatch.sourceforge.net/download/sleef-2.80.tar.gz"), sleef, SHA="d6c2314aee8c5201ca79237784aa2bb4fb3bbc727efa2f4386f04ad3d1bc0050")

prefix = usrdir(sleef)

#Hack for ordered dict -- don't create dict, just tuple of pairs
const simd_opts = (
"fma4" => ["-DENABLE_FMA4", "-mavx", "-mfma4"],
"avx2" => ["-DENABLE_AVX2", "-mavx2", "-mfma"],
"avx" => ["-DENABLE_AVX", "-mavx"],
"sse2" => ["-DENABLE_SSE2", "-msse2"],
# ARM NEON support is not quiet yet written
#"neon" => ["-DENABLE_NEON", "-mfloat-abi=softfp", "-mfpu=neon", "-static"]
)

const simd_flags = map(first, simd_opts)
const SIMD_ARCH = find_cpuflag(simd_flags, "purec", "JLSLEEF_ARCH")

if SIMD_ARCH=="purec"
	sleefsrcdir = joinpath(srcdir(sleef), "sleef-2.80/purec")
	const ARCHOPTS = String[]
	const prebuild = @build_steps begin
		`echo "building for No SIMD"`
	end
else
	sleefsrcdir = joinpath(srcdir(sleef), "sleef-2.80/simd")
	const ARCHOPTS = Dict(simd_opts)[SIMD_ARCH]
	const prebuild = @build_steps begin
		`echo "building for $SIMD_ARCH SIMD"`
		`ln -s sleefsimddp.c sleefdp.c`
		`ln -s sleefsimdsp.c sleefsp.c`
	end

end

const CCOPTS = split("-O -Wall -Wno-unused -Wno-attributes -lm")

if is_apple()
	const LIBBUILTOPTS = `-dynamiclib`
	const LIBFN = "sleef.dylib"
elseif is_linux()
	const LIBBUILTOPTS = `-shared`
	const LIBFN = "sleef.so"
else
	@assert is_windows()
	error("Windows not currently supported.")
end

provides(SimpleBuild,
    (@build_steps begin
        GetSources(sleef)
        @build_steps begin
            ChangeDirectory(sleefsrcdir)
            FileRule(joinpath(prefix,"lib", LIBFN),
                @build_steps begin
                    CreateDirectory(joinpath(prefix, "lib"))
					prebuild
					CCompile("sleefdp.c", "sleefdp.o", [CCOPTS; "-c"; "-fPIC"; ARCHOPTS],String[])
					CCompile("sleefsp.c", "sleefsp.o", [CCOPTS; "-c"; "-fPIC"; ARCHOPTS],String[])
					`gcc $LIBBUILTOPTS -o $LIBFN sleefsp.o sleefdp.o`
					#CCompile("sleefsimddp.o sleefsimdsp.o",LIBFN, [LIBBUILTOPTS]], String[] )
					`cp $LIBFN $(joinpath(prefix,"lib",LIBFN))`
                end
				)
        end
    end), sleef)

@BinDeps.install Dict([(:sleef, :jl_sleef)])

# Write the SIMD_ARCH to the deps.jl file.
# Recording what SIMD_ARCH was it build for is part of the deps
# and we need to make that available to the dep's consuming code
open(joinpath(BinDeps.depsdir(sleef), "deps.jl"), "a") do fp
	println(fp, "\nconst SIMD_ARCH = \"$(SIMD_ARCH)\"")
end


text:
# # Polarizability using automatic differentiation
#
# Simple example for computing properties using (forward-mode)
# automatic differentation.
# For a more classical approach and more details about computing polarizabilities,
# see [Polarizability by linear response](@ref).

using DFTK
using LinearAlgebra
using ForwardDiff

## Construct PlaneWaveBasis given a particular electric field strength
## Again we take the example of a Helium atom.
function make_basis(ε::T; a=10., Ecut=30) where T
    lattice=T(a) * I(3)  # lattice is a cube of ``a`` Bohrs
    ## Helium at the center of the box
    atoms     = [ElementPsp(:He, psp=load_psp("hgh/lda/He-q2"))]
    positions = [[1/2, 1/2, 1/2]]

    model = model_DFT(lattice, atoms, positions, [:lda_x, :lda_c_vwn];
                      extra_terms=[ExternalFromReal(r -> -ε * (r[1] - a/2))],
                      symmetries=false)
    PlaneWaveBasis(model; Ecut, kgrid=[1, 1, 1])  # No k-point sampling on isolated system
end

## dipole moment of a given density (assuming the current geometry)
function dipole(basis, ρ)
    @assert isdiag(basis.model.lattice)
    a  = basis.model.lattice[1, 1]
    rr = [a * (r[1] - 1/2) for r in r_vectors(basis)]
    sum(rr .* ρ) * basis.dvol
end

## Function to compute the dipole for a given field strength
function compute_dipole(ε; tol=1e-8, kwargs...)
    scfres = self_consistent_field(make_basis(ε; kwargs...), tol=tol)
    dipole(scfres.basis, scfres.ρ)
end;

# With this in place we can compute the polarizability from finite differences
# (just like in the previous example):
polarizability_fd = let
    ε = 0.01
    (compute_dipole(ε) - compute_dipole(0.0)) / ε
end

# We do the same thing using automatic differentiation. Under the hood this uses
# custom rules to implicitly differentiate through the self-consistent
# field fixed-point problem.
polarizability = ForwardDiff.derivative(compute_dipole, 0.0)
println()
println("Polarizability via ForwardDiff:       $polarizability")
println("Polarizability via finite difference: $polarizability_fd")


text:
module TestIsLess
    using Base.Test
    using CategoricalArrays

    pool = CategoricalPool([1, 2, 3])

    v1 = CategoricalValue(1, pool)
    v2 = CategoricalValue(2, pool)
    v3 = CategoricalValue(3, pool)

    @test_throws Exception v1 < v1
    @test_throws Exception v1 < v2
    @test_throws Exception v1 < v3
    @test_throws Exception v2 < v1
    @test_throws Exception v2 < v2
    @test_throws Exception v2 < v3
    @test_throws Exception v3 < v1
    @test_throws Exception v3 < v2
    @test_throws Exception v3 < v3

    @test_throws Exception v1 <= v1
    @test_throws Exception v1 <= v2
    @test_throws Exception v1 <= v3
    @test_throws Exception v2 <= v1
    @test_throws Exception v2 <= v2
    @test_throws Exception v2 <= v3
    @test_throws Exception v3 <= v1
    @test_throws Exception v3 <= v2
    @test_throws Exception v3 <= v3

    @test_throws Exception v1 > v1
    @test_throws Exception v1 > v2
    @test_throws Exception v1 > v3
    @test_throws Exception v2 > v1
    @test_throws Exception v2 > v2
    @test_throws Exception v2 > v3
    @test_throws Exception v3 > v1
    @test_throws Exception v3 > v2
    @test_throws Exception v3 > v3

    @test_throws Exception v1 >= v1
    @test_throws Exception v1 >= v2
    @test_throws Exception v1 >= v3
    @test_throws Exception v2 >= v1
    @test_throws Exception v2 >= v2
    @test_throws Exception v2 >= v3
    @test_throws Exception v3 >= v1
    @test_throws Exception v3 >= v2
    @test_throws Exception v3 >= v3

    @test ordered!(pool, true) === true
    @test ordered(pool) === true

    @test (v1 < v1) === false
    @test (v1 < v2) === true
    @test (v1 < v3) === true
    @test (v2 < v1) === false
    @test (v2 < v2) === false
    @test (v2 < v3) === true
    @test (v3 < v1) === false
    @test (v3 < v2) === false
    @test (v3 < v3) === false

    @test (v1 <= v1) === true
    @test (v1 <= v2) === true
    @test (v1 <= v3) === true
    @test (v2 <= v1) === false
    @test (v2 <= v2) === true
    @test (v2 <= v3) === true
    @test (v3 <= v1) === false
    @test (v3 <= v2) === false
    @test (v3 <= v3) === true

    @test (v1 > v1) === false
    @test (v1 > v2) === false
    @test (v1 > v3) === false
    @test (v2 > v1) === true
    @test (v2 > v2) === false
    @test (v2 > v3) === false
    @test (v3 > v1) === true
    @test (v3 > v2) === true
    @test (v3 > v3) === false

    @test (v1 >= v1) === true
    @test (v1 >= v2) === false
    @test (v1 >= v3) === false
    @test (v2 >= v1) === true
    @test (v2 >= v2) === true
    @test (v2 >= v3) === false
    @test (v3 >= v1) === true
    @test (v3 >= v2) === true
    @test (v3 >= v3) === true

    levels!(pool, [2, 3, 1])

    @test (v1 < v1) === false
    @test (v1 < v2) === false
    @test (v1 < v3) === false
    @test (v2 < v1) === true
    @test (v2 < v2) === false
    @test (v2 < v3) === true
    @test (v3 < v1) === true
    @test (v3 < v2) === false
    @test (v3 < v3) === false

    @test (v1 <= v1) === true
    @test (v1 <= v2) === false
    @test (v1 <= v3) === false
    @test (v2 <= v1) === true
    @test (v2 <= v2) === true
    @test (v2 <= v3) === true
    @test (v3 <= v1) === true
    @test (v3 <= v2) === false
    @test (v3 <= v3) === true

    @test (v1 > v1) === false
    @test (v1 > v2) === true
    @test (v1 > v3) === true
    @test (v2 > v1) === false
    @test (v2 > v2) === false
    @test (v2 > v3) === false
    @test (v3 > v1) === false
    @test (v3 > v2) === true
    @test (v3 > v3) === false

    @test (v1 >= v1) === true
    @test (v1 >= v2) === true
    @test (v1 >= v3) === true
    @test (v2 >= v1) === false
    @test (v2 >= v2) === true
    @test (v2 >= v3) === false
    @test (v3 >= v1) === false
    @test (v3 >= v2) === true
    @test (v3 >= v3) === true

    @test ordered!(pool, false) === false
    @test ordered(pool) === false

    @test_throws Exception v1 < v1
    @test_throws Exception v1 < v2
    @test_throws Exception v1 < v3
    @test_throws Exception v2 < v1
    @test_throws Exception v2 < v2
    @test_throws Exception v2 < v3
    @test_throws Exception v3 < v1
    @test_throws Exception v3 < v2
    @test_throws Exception v3 < v3

    @test_throws Exception v1 <= v1
    @test_throws Exception v1 <= v2
    @test_throws Exception v1 <= v3
    @test_throws Exception v2 <= v1
    @test_throws Exception v2 <= v2
    @test_throws Exception v2 <= v3
    @test_throws Exception v3 <= v1
    @test_throws Exception v3 <= v2
    @test_throws Exception v3 <= v3

    @test_throws Exception v1 > v1
    @test_throws Exception v1 > v2
    @test_throws Exception v1 > v3
    @test_throws Exception v2 > v1
    @test_throws Exception v2 > v2
    @test_throws Exception v2 > v3
    @test_throws Exception v3 > v1
    @test_throws Exception v3 > v2
    @test_throws Exception v3 > v3

    @test_throws Exception v1 >= v1
    @test_throws Exception v1 >= v2
    @test_throws Exception v1 >= v3
    @test_throws Exception v2 >= v1
    @test_throws Exception v2 >= v2
    @test_throws Exception v2 >= v3
    @test_throws Exception v3 >= v1
    @test_throws Exception v3 >= v2
    @test_throws Exception v3 >= v3
end


text:
export instruct, response, addbreak_every, show_cross, @read_args,
  @read_debug_args, oddball_paradigm, randomize_by
using ArgParse
using Juno: input, selector
import Juno

"""
    response(key1 => response1,key2 => response2,...;kwds...)

Create a watcher moment that records press of `key[n]` as
`record(response[n];kwds...)`.

See [`record`](@ref) for more details on how events are recorded.

When a key is pressed down, the record event occurs. Key releases are also
recorded, but are suffixed, by default, with "_up". This suffix can be changed
using the `keyup_suffix` keyword argument.
"""
function response(responses::Pair...;keyup_suffix="_up",info...)
  begin (event) ->
    for (key,response) in responses
      if iskeydown(event,key)
        record(response;info...)
      elseif iskeyup(event,key)
        record(response*keyup_suffix;info...)
      end
    end
  end
end

"""
    instruct(str;keys...)

Presents some instructions to the participant.

This adds "(Hit spacebar to continue...)" to the end of the text, and waits for
the participant to press spacebar to move on. It records an "instructions"
event to the data file.

Any keyword arguments are passed onto to `visual`, which can be used
to adjust how the instructions are displayed.
"""
function instruct(str;keys...)
  text = visual(str*" (Hit spacebar to continue...)";keys...)
  m = moment() do
    record("instructions")
    display(text)
  end
  [m,await_response(iskeydown(key":space:"))]
end

"""
    addbreak_every(n,total,
                   [response=key":space:"],[response_str="the spacebar"])

Adds a break every `n` times this event is added given a known number of
total such events.

By default this waits for the user to hit spacebar to move on.
"""
function addbreak_every(n,total,response=key":space:",
                        response_str="the spacebar")
  meta = Weber.metadata()
  index = meta[:break_every_index] = get(meta,:break_every_index,0) + 1
  if n <= index < total && (n == 1 || index % n == 1)
    message = moment() do
      display(visual("You can take a break. Hit "*
                     "$response_str when you're ready to resume... "*
                     "$(div(index,n)) of $(div(total-1,n)) breaks."))
    end

    addbreak(message,await_response(e -> iskeydown(e,response)))
  end
end

"""
    oddball_paradigm(trial_body_fn,n_oddballs,n_standards;
                     lead=20,n_standard_after_oddball=1,rng=Base.GLOBAL_RNG)

Helper to generate trials for an oddball paradigm.

The trial_body_fn should setup stimulus presentation: it takes one argument,
indicating if the stimulus should be a standard (false) or oddball (true)
stimulus.

It is usually best to use oddball_paradigm with a do block syntax. For instance,
the following code sets up 20 oddball and 150 standard trials.

    oddball_paradigm(20,150) do isoddball
      if isoddball
        addtrial(...create oddball trial here...)
      else
        addtrial(...create standard trial here...)
      end
    end

Alternatively, because `oddball_paradigm` returns the result of each
function call, you can pass the oddball paradigm as a series of moments
to a single trial.

    create_oddball_moment(isoddball) =
      isoddball ? moment(...oddball...) : moment(...standard...)
    addtrial(oddball_paradigm(create_oddball_moment,20,150))

# Keyword arguments

* **lead**: determines the number of standards that repeat before any oddballs
  get presented
* **oddball_spacing**: determines the number of standards after
  an oddball that must occur before a new oddball can occur
"""
function oddball_paradigm(fn,n_oddballs,n_standards;lead=20,oddball_spacing=1,
                          rng=Base.GLOBAL_RNG)
  min_standards = oddball_spacing * n_oddballs + lead
  @assert min_standards <= n_standards """
  You need at least $min_standards standards to have $n_oddballs oddballs
  because you need $oddball_spacing standards between each oddball and a lead of
  $lead standards.
  """

  oddballs_left = n_oddballs
  standards_left = n_standards
  n_stimuli = n_oddballs + n_standards
  n_last_standard = 0
  map(1:n_stimuli) do trial
    stimuli_left = n_stimuli - trial + 1
    min_trials_for_oddballs = oddballs_left*oddball_spacing
    oddball_chance = oddballs_left / (stimuli_left - min_trials_for_oddballs)

    if (trial > lead &&
        n_last_standard >= oddball_spacing &&
        rand(rng) < oddball_chance)

      n_last_standard = 0
      oddballs_left -= 1
      fn(true)
    else
      n_last_standard += 1
      standards_left -= 1
      fn(false)
    end
  end
end

"""
    show_cross([delta_t])

Creates a moment that shows a cross hair `delta_t` seconds after the start
of the previous moment (defaults to 0 seconds).
"""
function show_cross(delta_t::Number=0s;render_options...)
  moment(delta_t,display,"+";render_options...)
end

function as_arg(expr)
  if !isexpr(expr,:(=))
    error("Expected keyword parameters specifying additional program arguments.")
  end

  if expr.args[2] isa Symbol
    quote
      $(string(expr.args[1]))
      $(esc(:required)) = true
      $(esc(:arg_type)) = $(expr.args[2])
    end
  elseif expr.args[2].head == :vect
    quote
      $(string(expr.args[1]))
      $(esc(:required)) = true
      $(esc(:arg_type)) = String
      $(esc(:help)) = $(join(map(x -> x.args[1],expr.args[2].args),", "," or "))
    end
  else
    error("Expected keyword value to be a vector of symbols or a type.")
  end
end

function as_arg_checker(expr)
  if !(expr.args[2] isa Symbol) && expr.args[2].head == :vect
    quote
      let str = $(string(expr.args[1])), vals = $(expr.args[2])
        if !any(s -> string(s) == parsed[str],vals)
          println("Expected \"$str\" argument to be "*join(vals,", "," or ")*".")
          println(usage_string(s))
          exit()
        end
      end
    end
  else
    :nothing
  end
end

function as_arg_result(expr)
  if !(expr.args[2] isa Symbol) && expr.args[2].head == :vect
    :(Symbol(parsed[$(string(expr.args[1]))]))
  else
    :(parsed[$(string(expr.args[1]))])
  end
end

"""
    randomize_by(itr)

Randomize by a given iterable object, usually a string (e.g. the subject id.)

If the same iterable is given, calls to random functions (e.g. `rand`, `randn`
and `shuffle`) will result in the same output.
"""
randomize_by(itr) = srand(reinterpret(UInt32,collect(itr)))

"""
    @read_debug_args(description,[keyword args...])

Same as @read_args, but better suited to debugging errors in your program
when running the experiment in Juno.

Specifically, this verison will never spawn a new process to run the
experiment. This means that you can safely step through the code using debugging
tools. In this case, you will also likely want to set `debug=true` when defining
your [`Experiment`](@ref) object.
"""
macro read_debug_args(description,keys...)
  expr = :(@read_args($description))
  for key in keys
    push!(expr.args,key)
  end
  quote
    Weber.start_debug_read_args()
    result = $expr
    Weber.end_debug_read_args()
    result
  end
end

"""
    @read_args(description,[keyword args...])

Reads experimental parameters from the user.

With no additional keyword arguments this requests the subject id, and an
optional `skip` parameter (defaults to 0) from the user, and then returns them
both in a tuple. The skip can be used to restart an experiment by passing it as
the `skip` keyword argument to the `Experiment` constructor.

You can specify additional keyword arguments to request additional
values from the user. Arguments that are a type will yield a request for
textual input, and will verify that that input can be parsed as the given type.
Arguments whose values are a list of symbols yield a request that the user select
one of the specified values.

Arguments are requested from the user either as command-line arguments,
or, if no command-line arguments were specified, interactively. Interactive
arguments work both in the terminal or in Juno. This macro also
generates useful help text that will be displayed to the user
when they give a single command-line "-h" argument. This help text
will include the `desecription` string.

# Example

    subject_id,skip,condition,block = @read_args("A simple experiment",
      condition=[:red,:green,:blue],block=Int)
"""
macro read_args(description,keys...)
  arg_expr = quote
    "sid"
    $(esc(:help)) = "The subject id."
    $(esc(:required)) = true
    $(esc(:arg_type)) = String
  end

  for arg_body in map(as_arg,keys)
    for line in arg_body.args
      push!(arg_expr.args,line)
    end
  end

  skip_expr = quote
    "skip"
    $(esc(:help)) = "# of offsets to skip. Useful for restarting in middle of experiment."
    $(esc(:required)) = false
    $(esc(:arg_type)) = Int
    $(esc(:default)) = 0
  end
  arg_expr.args = vcat(arg_expr.args,skip_expr.args)

  arg_body = quote
    s = ArgParseSettings(description = $(esc(description)))

    @add_arg_table s begin
      $arg_expr
    end

    parsed = parse_args(ARGS,s)
  end

  for line in map(as_arg_checker,keys)
    push!(arg_body.args,line)
  end

  result_tuple = :((parsed["sid"],parsed["skip"]))
  for result = map(as_arg_result,keys)
    push!(result_tuple.args,result)
  end
  push!(arg_body.args,result_tuple)

  script_file = gensym(:script_file)
  collect_args = :(collect_args($(esc(description)),$script_file))
  for key in keys
    if !isexpr(key,:(=))
      error("Expected keyword parameters specifying additional program arguments.")
    else
      push!(collect_args.args,Expr(:kw,key.args...))
    end
  end

  quote
    cd(dirname(@__FILE__))
    if length(ARGS) > 0
      $arg_body
    else
      $script_file = @__FILE__
      $collect_args
    end
  end
end

start_debug_read_args() = global debug_read_args = true
end_debug_read_args() = global debug_read_args = false
debug_read_args = false
function collect_args(description,script_file;keys...)
  if Juno.isactive()
    info("Please type your responses directly below the prompt."*
         "Entering them further down (next to the '>') will not work.")
  end

  print("Enter subject id: ")
  sid = chomp(input())
  args = Array{Any}(length(keys))
  for (i,(kw,value)) in enumerate(keys)
    if value isa Type
      except = true
      while except
        try
          print("Enter $kw: ")
          args[i] = parse(value,chomp(input()))
          except = false
        catch
          except = true
          println("Expected $value.")
        end
      end
    else
      error = true
      while error
        print("Enter $kw ($(join(map(string,value),", "," or "))): ")
        args[i] = chomp(input())
        if Symbol(args[i]) ∉ value
          println("Expected $kw to be $(join(map(string,value),", "," or ")) "*
                  "but got $(args[i]).")
          error = true
        else
          error = false
        end
      end
    end
  end
  print("Offset to start at? (default = 0): ")
  str = input()
  if isempty(chomp(str))
    skip = 0
  else
    skip = parse(Int,str)
  end

  if Juno.isactive() && debug_read_args == false
    info("Spawning experiment in new child process...")
    info("If you need to debug your experiment, use @read_debug_args instead.")
    timestr = Dates.format(now(),"yyyy-mm-dd__HH_MM_SS")
    logfile = "weber_$timestr.log"
    child = spawn(pipeline(`$(joinpath(JULIA_HOME,"julia")) $script_file $sid $args $skip`,
                           stdout=logfile,stderr=logfile))
    info("Experiment log will be written to $logfile. The experiment should begin shortly...")
    wait(child)
    info("Experiment complete. Terminating julia, so the next run of an "*
         "expeirment runs cleanly.")
    exit()
  else
    info("Running...")
    (sid,skip,args...)
  end
end


text:
@testset "no intercept model" begin
    x = [-5.0, -4.8, -4.6, -4.4, -4.2, -4.0, -3.8, -3.6, -3.4, -3.2, -3.0, -2.8, -2.6, -2.4, -2.2, -2.0, -1.8, -1.6, -1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0]
    y = [32.76095293186829, 32.114646881731815, 30.482249790567508, 30.172132123809526, 28.858546806301607, 27.55525364876047, 26.748316399220183, 25.47632057519211, 24.178757213602406, 23.115113146604088, 21.67768685139419, 20.667430892918777, 19.25030326250073, 17.969338480523586, 16.92778066296735, 16.01093498430364, 14.972217418295738, 14.129273625410654, 12.528548495622791, 11.995655292690893, 10.437691863433523, 9.04875056404323, 8.352545547062054, 6.595027421481294, 5.769353744384407, 4.656766397706274, 3.2769464299313293, 2.557398138612538, 1.8697655690559674, 0.1758423094007262, -1.1812180900645513, -2.4595328718334226, -3.224902727206241, -4.52850811714283, -5.620262454732497, -6.65705617649476, -8.109350054557144, -8.964028166630365, -9.764063784932889, -11.49260887582183, -12.144819717782108, -13.843175083307838, -14.603379688823846, -15.335111221096993, -16.28627855774275, -17.73382957392828, -19.156926371425502, -20.641255379884615, -21.17362575264, -22.5708673184868, -23.072613984429083]
    df = DataFrame(x=x, y=y)
    lrnoint = regress(@formula(y ~ 0 + x), df, req_stats=[:default, :aic, :vif])
    @test isapprox(0.9252939081, lrnoint.R2)
    @test isapprox(0.9237997863, lrnoint.ADJR2)
    @test isapprox(160.90768623, lrnoint.AIC)
    @test isapprox([0.], lrnoint.VIF)

    results = predict_in_sample(lrnoint, df, α=0.05, req_stats=["default", "rstudent"])
    @test isapprox([28.385804339677044, 27.25037216608996, 26.114939992502876, 24.9795078189158, 23.84407564532872, 22.708643471741635, 21.57321129815455, 20.43777912456747, 19.30234695098039, 18.16691477739331, 17.031482603806225, 15.896050430219143, 14.760618256632062, 13.62518608304498, 12.4897539094579, 11.354321735870817, 10.218889562283735, 9.083457388696655, 7.948025215109571, 6.81259304152249, 5.677160867935409, 4.541728694348327, 3.406296520761245, 2.2708643471741636, 1.1354321735870818, 0.0, -1.1354321735870818, -2.2708643471741636, -3.406296520761245, -4.541728694348327, -5.677160867935409, -6.81259304152249, -7.948025215109571, -9.083457388696655, -10.218889562283735, -11.354321735870817, -12.4897539094579, -13.62518608304498, -14.760618256632062, -15.896050430219143, -17.031482603806225, -18.16691477739331, -19.30234695098039, -20.43777912456747, -21.57321129815455, -22.708643471741635, -23.84407564532872, -24.9795078189158, -26.114939992502876, -27.25037216608996, -28.385804339677044]
            , results.predicted)
    @test isapprox([0.938033258871313, 1.0426179367896207, 0.9319648731276958, 1.1097378448312425, 1.0685340559242513, 1.0299419326232258, 1.0993609703148552, 1.0677740527805426, 1.0309373482675432, 1.0448144116724511, 0.978326230602586, 1.0038487427946268, 0.9422793707140289, 0.9101464457061582, 0.9291611309816967, 0.9748389070719975, 0.9946108960924838, 1.0563303532318353, 0.9563424123040439, 1.0843302339300107, 0.993564396349128, 0.9392838550370636, 1.032388388246509, 0.9000336054567951, 0.9655586062937045, 0.97036522335683, 0.9185810103197727, 1.007013102837162, 1.1028816259631355, 0.9840106274248026, 0.9373227029099432, 0.9074862070811145, 0.9866973636129569, 0.9515650780179017, 0.9616191528659762, 0.9835165279774034, 0.9168886277207625, 0.9778135787709513, 1.0509330791271387, 0.9250288747899885, 1.0300270784500691, 0.910564302144073, 0.992654638302146, 1.0816862270131862, 1.1237223658970572, 1.0578001520839455, 0.9973186697257519, 0.9236477076314706, 1.0570827343881202, 1.0021858611232244, 1.1440369678715858]
            , results.rstudent)

    y = [3.547744106900422, 9.972950249405148, 16.471345464154027, 22.46768807351274, 20.369933318011807, 21.18590757820348, 29.962620198209024, 30.684400502954748, 29.28429078492597, 34.272759386588824, 33.05504692986838, 45.09273876302829, 45.28374262744938, 54.54563960566191, 46.86173948296966, 46.85926120310666, 67.54337216713414, 66.0400205145086, 64.77001443647681, 63.98759256558095, 15.016939388490687, 12.380885701920008, 10.221963402745288, 24.826987790646672, 22.231511892187548, 18.05125492502642, 21.809661866284717, 28.533306224702308, 24.70476800009685, 39.592181057942916, 37.425282624708906, 33.89372811236659, 44.07442335927311, 42.61943719318272, 52.32531447897055, 44.12096163880534, 49.19965880584543, 52.304114239221036, 59.79488104919937, 65.00419916894333, 10.978506745605673, 12.944637189081874, 16.96525561080181, 15.93489355798633, 33.559135307088305, 28.887571687420433, 28.296418324622593, 32.2405215537489, 31.466024917490223, 37.78855308849255, 48.09725215994402, 40.70190542438069, 43.2525436240948, 44.75159242105558, 52.553066965996074, 53.265851121437095, 60.092700204726015, 62.78014347092756, 69.07895282754228, 74.27890668517966, 10.57283140105129, 7.290600131361426, 16.10299402050687, 14.428954773831242, 22.418180226673464, 28.21022852582732, 29.60436622143203, 33.648588929669636, 37.451930576147994, 43.85548900583812, 41.44168340404242, 43.48815266671309, 54.72835160956764, 53.55177468229062, 50.62088969800169, 50.863408563713335, 47.40251347184729, 64.29413401802115, 64.20687412126479, 73.66653216085827, 15.820723594639162, 21.973463928234743, 23.804440715385162, 15.139822408433913, 30.015089890369985, 28.08454394421385, 33.041880065463566, 28.49429418531324, 38.33763660905526, 34.503176013521724, 48.748946870573235, 45.45351516824085, 53.522908096188765, 45.95216131836423, 63.13018379633756, 63.4236036208151, 65.28265579397677, 55.43500146374922, 76.50187470137375, 67.22421998359037, 19.289152114521762, 27.63557010588222, 17.700031078096686, 25.462368278660957, 22.94613580060685, 36.19731649621813, 35.22216995579936, 25.526434727578838, 45.882864925557726, 38.71433797181679, 50.617762276434554, 41.96951039650285, 52.265123328453214, 45.383991138243765, 62.7923270665056, 64.40670612696276, 74.89775274405821, 72.89056716347118, 66.37071343209973, 76.32913721410918, 29.923781174452596, 16.465832617450324, 30.530915733275403, 30.512411156491954, 28.04012351007911, 26.315140074869376, 37.18491928428231, 41.958551085353626, 48.370736387628895, 49.419917216561835, 48.67296268029715, 55.484477112881166, 51.120639229597025, 54.797092949987224, 61.92608065418044, 69.79495109420618, 64.43521939892794, 74.35280205157255, 77.22355341160723, 68.90654715174155, 10.87752158238755, 20.25006748725279, 31.71957876614495, 31.42898857400639, 32.505052996368256, 34.1225641368903, 35.66496173153403, 31.76371648643483, 42.97107454773545, 47.41827763710622]
    x2 = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8]
    x3 = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    df = DataFrame(y=y, x2=x2, x3=x3)

    lrnoint = regress(@formula(y ~ 0 + x2 +x3), df, req_stats=["default", "pcorr1", "pcorr2", "scorr1", "scorr2", "vif"])
    @test isapprox([207017.27329467665, 84364.93565847247], lrnoint.Type1SS)
    @test isapprox([6502.255171530228, 84364.93565847244], lrnoint.Type2SS)
    @test isapprox([0.9814816505541067, 0.9557504161506739], lrnoint.pcorr1)
    @test isapprox([0.6247239646866073, 0.9557504161506739], lrnoint.pcorr2)
    @test isapprox([0.7010686580206047, 0.28570375449728647], lrnoint.scorr1)
    @test isapprox([0.02202003356851988, 0.28570375449728636], lrnoint.scorr2)
    @test isapprox([1.0099808403511108, 1.0099808403511117], lrnoint.VIF)
    @test isapprox(0.9867724125178912, lrnoint.R2)
    @test isapprox(0.9865936613357005, lrnoint.ADJR2)
end

@testset "no intercept weighted regression" begin
        tw = [
            2.3  7.4  0.058
            3.0  7.6  0.073
            2.9  8.2  0.114
            4.8  9.0  0.144
            1.3 10.4  0.151
            3.6 11.7  0.119
            2.3 11.7  0.119
            4.6 11.8  0.114
            3.0 12.4  0.073
            5.4 12.9  0.035
            6.4 14.0  0
        ] # data from https://blogs.sas.com/content/iml/2016/10/05/weighted-regression.html

        df = DataFrame(tw, [:y,:x,:w])
        f = @formula(y ~ 0 + x)
        lm = regress(f, df, weights="w")

        @test isapprox([0.3056575637534476], lm.coefs)
        @test isapprox(0.8626294380695502, lm.R2)
        @test isapprox(0.8473660422995002, lm.ADJR2)
        @test isapprox([0.040658241629241754], lm.stderrors)
        @test isapprox([3.6247531175150984e-5], lm.p_values)

end


text:
import Base:get, length, show

export CosDict, CosString, CosNumeric, CosBoolean, CosTrue, CosFalse,
       CosObject, CosNull, CosNullType,CosFloat, CosInt, CosArray, CosName,
       CosDict, CosIndirectObjectRef, CosStream, set!, @cn_str,
    createTreeNode, CosTreeNode, CosIndirectObject,
    CosDictType

"""
```
    CosObject
```
PDF is a structured document format with lots of internal data structures like
dictionaries, arrays, trees. `CosObject` is the interface to access these objects and get
detailed access to the objects and gather additional information. Although, defined in the
COS layer, objects of these type are returned from almost all the APIs. Hence, the objects
have a separate significance whether you need to use the `Cos` layer or not. Below is the
object hierarchy.

```
CosObject                           Abstract
    CosNull                         Value (CosNullType)
CosString                           Abstract
CosName                             Concrete
CosNumeric                          Abstract
    CosInt                          Concrete
    CosFloat                        Concrete
CosBoolean                          Concrete
    CosTrue                         Value (CosBoolean)
    CosFalse                        Value (CosBoolean)
CosDict                             Concrete
CosArray                            Concrete
CosStream                           Concrete (always wrapped as an indirect object)
CosIndirectObjectRef                Concrete (only useful when CosDoc is available)
```
"""
abstract type CosObject end

get(o::T) where {T <: CosObject} = o.val

"""
```
    CosString
```
Abstract type that represents a PDF string. In PDF objects are mere byte representations.
They translate to actual text strings by application of fonts and associated encodings.
"""
abstract type CosString <: CosObject end

get(o::T) where {T <: CosString} = copy(o.val)
"""
```
    CosNumeric
```
Abstract type for numeric objects. The objects can be an integer [`CosInt`](@ref) or float
[`CosFloat`](@ref).
"""
abstract type CosNumeric <: CosObject end

"""
```
    CosBoolean
```
A boolean object in PDF which is either a `CosTrue` or `CosFalse`
"""
struct CosBoolean <: CosObject
    val::Bool
end

const CosTrue=CosBoolean(true)
const CosFalse=CosBoolean(false)

struct CosNullType <: CosObject end

"""
```
    CosNull
```
PDF representation of a `null` object. Can be applied to [`CosObject`](@ref)
of any type.
"""
const CosNull=CosNullType()

"""
```
    CosFloat
```
A numeric float data type.
"""
struct CosFloat <: CosNumeric
    val::Float32
end

"""
```
    CosInt
```
An integer in PDF document.
"""
struct CosInt <: CosNumeric
    val::Int
end

"""
```
    CosIndirectObjectRef
```
A parsed data structure to ensure the object information is stored as an object.
This has no meaning without a associated CosDoc. When a reference object is hit
the object should be searched from the CosDoc and returned.
"""
struct CosIndirectObjectRef <: CosObject
  val::Tuple{Int,Int}
  CosIndirectObjectRef(num::Int, gen::Int)=new((num,gen))
end

mutable struct CosIndirectObject{T <: CosObject} <: CosObject
  num::Int
  gen::Int
  obj::T
end

get(o::CosIndirectObject) = get(o.obj)

"""
```
    CosName
```
Name objects are symbols used in PDF documents.
"""
struct CosName <: CosObject
    val::Symbol
    CosName(str::AbstractString) = new(Symbol("CosName_",str))
end

"""
```
    @cn_str(str) -> CosName
```
A string decorator for easier instantiation of a [`CosName`](@ref)
"""
macro cn_str(str)
    return CosName(str)
end

"""
```
    CosXString
```
Concrete representation of a [`CosString`](@ref) object. The underlying data is
represented as hexadecimal characters in ASCII.
"""
struct CosXString <: CosString
  val::Vector{UInt8}
  CosXString(arr::Vector{UInt8})=new(arr)
end

"""
```
    CosLiteralString
```
Concrete representation of a [`CosString`](@ref) object. The underlying data is
represented by byte representations without any encoding.
"""
struct CosLiteralString <: CosString
    val::Vector{UInt8}
    CosLiteralString(arr::Vector{UInt8}) = new(arr)
end

function CosLiteralString(str::AbstractString)
    buf = IOBuffer()
    for c in str
        print(buf, Char(c))
    end
    CosLiteralString(take!(buf))
end

"""
```
    CosArray
```
An array in a PDF file. The objects can be any combination of [`CosObject`](@ref).
"""
mutable struct CosArray <: CosObject
    val::Vector{CosObject}
    function CosArray(arr::Vector{T} where {T<:CosObject})
      val = Vector{CosObject}()
      for v in arr
        push!(val,v)
      end
      new(val)
    end
    CosArray()=new(Vector{CosObject}())
end

"""
```
    get(o::CosArray, isNative=false) -> Vector{CosObject}
```
An array in a PDF file. The objects can be any combination of
[`CosObject`](@ref).

`isNative = true` will return the underlying native object inside the `CosArray`
by invoking get method on it.
"""
get(o::CosArray, isNative=false) = isNative ? map((x)->get(x), o.val) : o.val

get(o::CosIndirectObject{CosArray}, isNative=false) = get(o.obj, isNative)
"""
```
    length(o::CosArray) -> Int
```
Length of the `CosArray`
"""
length(o::CosArray) = length(o.val)

length(o::CosIndirectObject{CosArray}) = length(o.obj)

"""
```
    CosDict
```
Name value pair of a PDF objects. The object is very similar to the `Dict`
object. The `key` has to be of a [`CosName`](@ref) type.
"""
mutable struct CosDict <: CosObject
    val::Dict{CosName, CosObject}
    CosDict()=new(Dict{CosName, CosObject}())
end

const CosDictType = Union{CosDict, CosIndirectObject{CosDict}}

"""
```
    get(dict::CosDict, name::CosName) -> CosObject
```
Returns the value as a [`CosObject`](@ref) for the key `name`
"""
get(dict::CosDict, name::CosName) = get(dict.val, name, CosNull)

get(o::CosIndirectObject{CosDict}, name::CosName) = get(o.obj, name)

"""
```
    set!(dict::CosDict, name::CosName, obj::CosObject) -> CosObject
```
Sets the value on a dictionary object. Setting a `CosNull` object deletes the
object from the dictionary.
"""
function set!(dict::CosDict, name::CosName, obj::CosObject)
    if (obj === CosNull)
        delete!(dict.val,name)
    else
        dict.val[name] = obj
    end
    return dict
end

set!(o::CosIndirectObject{CosDict}, name::CosName, obj::CosObject) =
    set!(o.obj, name, obj)

"""
```
    CosStream
```
A stream object in a PDF. Stream objects have an `extends` disctionary, followed
by binary data.
"""
mutable struct CosStream <: CosObject
    extent::CosDict
    isInternal::Bool
    CosStream(d::CosDict,isInternal::Bool=true) = new(d, isInternal)
end

get(stm::CosStream, name::CosName) = get(stm.extent, name)

get(o::CosIndirectObject{CosStream}, name::CosName) = get(o.obj, name)

set!(stm::CosStream, name::CosName, obj::CosObject)=
    set!(stm.extent, name, obj)

set!(o::CosIndirectObject{CosStream}, name::CosName, obj::CosObject) =
    set!(o.obj, name, obj)

"""
Decodes the stream and provides output as an BufferedInputStream.
"""
get(stm::CosStream) = decode(stm)

"""
```
    CosObjectStream
```
"""
mutable struct CosObjectStream <: CosObject
    stm::CosStream
    n::Int
    first::Int
    oids::Vector{Int}
    oloc::Vector{Int}
    function CosObjectStream(s::CosStream)
        n = get(s, CosName("N"))
        @assert n != CosNull
        first = get(s, CosName("First"))
        @assert first != CosNull
        cosStreamRemoveFilters(s)
        n_n = get(n)
        first_n = get(first)
        oids = zeros(Int, n_n)
        oloc = zeros(Int, n_n)
        read_object_info_from_stm(s, oids, oloc, n_n, first_n)
        new(s, n_n, first_n,oids, oloc)
    end
end

get(os::CosObjectStream, name::CosName) = get(os.stm, name)

get(os::CosIndirectObject{CosObjectStream}, name::CosName) = get(os.obj,name)

set!(os::CosObjectStream, name::CosName, obj::CosObject)=
    set!(os.stm, name, obj)

set!(os::CosIndirectObject{CosObjectStream}, name::CosName, obj::CosObject)=
    set!(os.obj,name,obj)

get(os::CosObjectStream) = get(os.stm)

"""
```
    CosObjectStream
```
"""
mutable struct CosXRefStream<: CosObject
  stm::CosStream
  isDecoded::Bool
  function CosXRefStream(s::CosStream, isDecoded::Bool=false)
      new(s, isDecoded)
  end
end

get(os::CosXRefStream, name::CosName) = get(os.stm, name)

get(os::CosIndirectObject{CosXRefStream}, name::CosName) = get(os.obj,name)

set!(os::CosXRefStream, name::CosName, obj::CosObject)=
    set!(os.stm, name, obj)

set!(os::CosIndirectObject{CosXRefStream}, name::CosName, obj::CosObject)=
    set!(os.obj,name,obj)

get(os::CosXRefStream) = get(os.stm)

"""
Can be a Number Tree or a Name Tree.

`kids`: is `null` in case of a leaf node
`range`: is `null` in case of a root node
`values`: is `null` in case of an intermediate node

Intent: faster loookup without needing to load the complete tree structure.
Hence, the tree will not be loaded on full scan.
"""
mutable struct CosTreeNode{K <: Union{Int, String}}
    values::Union{Nothing, Vector{Tuple{K, CosObject}}}
    kids::Union{Nothing, Vector{CosIndirectObjectRef}}
    range::Union{Nothing, Tuple{K, K}}
    function CosTreeNode{K}() where {K <: Union{Int, String}}
        new(nothing, nothing, nothing)
    end
end

# If K is Int, it's a number tree else it's a String which is a name tree
function createTreeNode(::Type{K}, dict::CosObject) where K
    range = get(dict, CosName("Limits"))
    kids  = get(dict, CosName("Kids"))
    node = CosTreeNode{K}()
    if (range !== CosNull)
        r = get(range, true)
        node.range = (r[1], r[2])
    end
    if (kids !== CosNull)
        node.kids = get(kids)
    end
    return populate_values(node, dict)
end

function populate_values(node::CosTreeNode{Int}, dict::CosObject)
    nums = get(dict, CosName("Nums"))
    if (nums !== CosNull)
        v = get(nums)
        values = [(get(v[2i-1]), v[2i]) for i=1:div(length(v), 2)]
        node.values = values
    end
    return node
end

function populate_values(node::CosTreeNode{String}, dict::CosObject)
    names = get(dict, CosName("Names"))
    if (names !== CosNull)
        v = get(names, true)
        values = [(get(v[2i-1]), v[2i]) for i=1:div(length(v), 2)]
        node.values = values
    end
    return node
end

# All show methods

show(io::IO, o::CosObject) = print(io, o.val)

showref(io::IO, o::CosObject) = show(io, o)

show(io::IO, o::CosNullType) = print(io, "null")

show(io::IO, o::CosName) = print(io, "/", String(o))

show(io::IO, o::CosXString) =  print(io, "<", String(copy(o.val)), ">")

show(io::IO, o::CosLiteralString) = print(io, "(", String(copy(o.val)), ")")

function show(io::IO, o::CosArray)
  print(io, '[')
  for obj in o.val
    showref(io, obj)
    print(io, ' ')
  end
  print(io, ']')
end

function show(io::IO, o::CosDict)
    print(io, "<<\n")
    for (k,v) in o.val
        print(io, '\t')
        show(io, k)
        print(io, '\t')
        showref(io, v)
        println(io, "")
    end
    print(io, ">>")
end

show(io::IO, stm::CosStream) =
  (show(io, stm.extent); print(io, "\nstream\n...\nendstream"))

show(io::IO, os::CosObjectStream) = show(io, os.stm)

show(io::IO, o::CosIndirectObjectRef) = print(io, o.val[1], ' ', o.val[2], " R")

showref(io::IO, o::CosIndirectObject) = print(io, o.num, ' ', o.gen, " R")

function show(io::IO, o::CosIndirectObject)
    println(io, "")
    println(io, o.num, ' ', o.gen, " obj")
    print(io, o.obj)
    println(io, "\nendobj\n")
end

"""
```
    CosComment
```
A comment object in PDF which is normally ignored as a whitespace.
"""
struct CosComment <: CosObject
    val::String
end

CosComment(barr::Vector{UInt8}) = CosComment(String(Char.(barr)))

show(io::IO, os::CosComment) = print(io, '%', os.val)


text:
module Day01
export day01

function readInput(inputfile="inputs/input01.txt")
    lines = parse.(Int, eachline(inputfile))
    return lines
end

function day01(inputdata = readInput("inputs/input01.txt"))
    return (day01Part1(inputdata), day01Part2(inputdata))
end

function day01Part1(values)
    #count(isless.(values[1:end-1], values[2:end]))
    #@views mapreduce(isless, +, values[1:end-1], values[2:end])
    return lagged_compare(values, 1)
end

function day01Part2(values)
    #count(isless.(values[1:end-3], values[4:end]))
    #@views mapreduce(isless, +, values[1:end-3], values[4:end])
    return lagged_compare(values, 3)
end

function lagged_compare(values, lag=1)
    count = 0
    for i in 1:length(values)-lag
        if values[i+lag] > values[i]
            count += 1
        end
    end
    return count
end

end #module


text:
include(joinpath(@__DIR__, "es", "standard_es.jl"))

#retuns a vector of ints, each of that
function split_numbers_10³(num::Integer)
    digits(num,base=10^3)
end

function split_numbers_10⁶(num::Integer)
    digits(num,base=10^6)
end

#function
function es_spell_1e3(number,short_one=false)
    if number <= 99 #0 - 99
        if number<30 #1 - 30
            iszero(number) && return ""
            if isone(number) && short_one
                return "uno"
            elseif number == 21 && short_one
                return "veintiuno"
            else
                return es_1_a_29[number]
            end
        else #30-99
            unit,dec = digits(number)
            if isone(unit) && short_one
                unit_text = " y uno"
            elseif iszero(unit)
                unit_text = ""
            else
                unit_text =  " y " * es_1_a_29[unit]
            end
            return es_decenas[dec] * unit_text
        end
    elseif 100  <= number <= 999
        number == 100 && return "cien"
        unit,cent = digits(number,base=100)
        unit_text = es_spell_1e3(unit,short_one)
        return es_centenas[cent] * " "  * unit_text
    end
end

function es_spell_1e6(number,short_one=false)
    number = Int(number)
    number == 0 && return ""
    number == 1000 && return "mil"
    number <= 999 && return es_spell_1e3(number,short_one)
    low,high = digits(number,base=1000)
    low_txt = es_spell_1e3(low,short_one)
    if high == 1
        high_text = "mil "
    else
        high_text = es_spell_1e3(high,false) * " mil "
    end
    return high_text * low_txt
end

#spells only decimals, including 0
function es_spell_decimal(number,partitive=true)
    whole, decimal = split(string(number), ".")
    decnum,wholenum = abs.(modf(number))
    decimal_length = length(decimal)
    decint = decnum*10^decimal_length
    if partitive
        if decimal_length>6
            throw(error("""partitive pronunciation of spanish decimals is only supported to up to six decimals, the decimal part is $decimal .
            If you want english-like decimals (cero punto dos), try using partitive = false.
            If the decimal part seems to differ from the perceived input number, try using big numbers: `big"123.23"`.
            """))
        end

        dectext = es_spell_1e6(Dec128(decint),false)
        res = strip(dectext) * " " * es_decimales[decimal_length]
        if decint != 1
            res = res * "s"
        end
        return res
    else
        res =  digits(decint) .|> z->es_spell_1e3(z,true) .|> strip |> reverse |> z->join(z," ")
        res = "punto " * res
    end
end

function es_spell_large_map(number,i)
    if isone(i)
        return es_spell_1e6(number,true)
    elseif isone(number)
        return "un " * es_multiplos_1e6_singular[i-1]
    else

        return es_spell_1e6(number,false) * " " * es_multiplos_1e6_plural[i-1]
    end
end

function es_spell_large(_number)
    number = Int(abs(_number)) #could cause problems, as we cannot convert (yet) to big int
    list =digits(number,base=1_000_000)
    res = es_spell_large_map.(list,1:length(list)) .|> strip |> reverse |> z->join(z," ")
    (_number < 0) && (res = "menos " * res)
    return res
end

function spelled_out_es(number,dict=:standard)
    if iszero(number)
        return "cero"
    end
    if isinteger(number)
        if dict in (:standard,:large,:modern)
            res = es_spell_large(number)
        elseif dict in (:short,)
            res = es_spell_short(number)
        else
            throw(error("unrecognized dict value: $dict"))
        end
        return res
    else #with decimals

        #delete this when proper separate keywords are implemented
        partitive = true
        _dec,_int = modf(number)
        _int = Dec128(_int)
        if (_int == 0) & !partitive
            intres = "cero"

        elseif (_int == 0) & partitive
            intres = ""
        elseif dict in (:standard,:large,:modern)
            intres = es_spell_large(_int) * " con "
        elseif dict in (:short,)
            intres = es_spell_short(_int) * " con "
        else
            throw(error("unrecognized dict value: $dict"))
        end

        decres = es_spell_decimal(abs(_dec))
        return intres * decres
    end
end


text:
module DirichletFESpacesTests

using Test
using Gridap.Algebra
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.CellData

n = 2
mesh = (n,n)
domain = (0,1,0,1)
order = 1
dirichlet_tags = [1,3,4,6,7]
model = CartesianDiscreteModel(domain, mesh)
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"dirichlet",dirichlet_tags)

Vf = TestFESpace(
  reffe=:Lagrangian, order=order, valuetype=Float64,
  conformity=:H1, model=model, labels=labels, dirichlet_tags="dirichlet")

Vd = DirichletFESpace(Vf)
matvecdata = ([],[],[])
matdata = ([],[],[])
vecdata = ([],[])
test_single_field_fe_space(Vd,matvecdata,matdata,vecdata)

Uf = TrialFESpace(Vf,1)

Ud = TrialFESpace(Vd)
test_single_field_fe_space(Ud,matvecdata,matdata,vecdata)

trian = Triangulation(model)
degree = 2*order
quad = CellQuadrature(trian,degree)

a(u,v) = u*v
l(v) = v*0
t_Ω = AffineFETerm(a,l,trian,quad)

op_ff = AffineFEOperator(Uf,Vf,t_Ω)
op_fd = AffineFEOperator(Ud,Vf,t_Ω)
op_df = AffineFEOperator(Uf,Vd,t_Ω)

xd = get_dirichlet_values(Uf)

@test maximum(abs.(get_vector(op_ff) + get_matrix(op_fd)*xd)) < 1.0e-10

end # module


text:
# GENERATOR PENALTY MODEL FORMULATION THAT PENALIZES VIOLATION OF GENERATOR BOUNDS INSTEAD OF USING RESERVE SATURATION

using Distributions, StatsFuns, StatsBase, JuMP, Gurobi, Ipopt

#*========= MODEL INFORMATION ==========
const model = "118bus"
const modelFile = model * ".jl"
include(modelFile)
#*======================================

include("utilities.jl")
include("modelFns.jl")
include("solveRecourse.jl")
include("evaluateSoln.jl")
include("projection.jl")

#*============ GENERATOR BOUNDS PENALTY ==========

# penalty parameters for generator bounds violation
const gamma_gen_val = logspace(0.0,5.0,16)
if(model == "118bus")
	gamma_gen_val = logspace(0.0,4.0,9)
end
const delta_gen = 1.0

#*============ LINE PENALTY ==========
# set line penalty coefficients to get comparable results
const gamma_line_val = logspace(1.0,5.0,9)

#*============ REPLICATES AND SCENARIOS ==========

# number of replicates for each parameter setting
const numReplicates = 5
const startingRep = 1

if(model == "118bus")
	numReplicates = 1
end

# number of scenarios for the SAA problems
const numScenarios = 500

# set seed for reproducibility
srand(4428)

# generate and store Gurobi environment once and for all
const gurobi_env = Gurobi.Env()

# solve the GP SAA model (convex QP) with penalties on generator bounds violation
function solveGPModel(gamma_gen::Float64,delta_gen::Float64,gamma_line::Float64,delta_line::Float64,
				omega_d::Array{Float64},omega_w::Array{Float64},useGurobi::Bool=true)

	# compute sum of demand deviations for each scenario
	sumMeanDemands = sum(meanDemand[j] for j=1:numLoads)
	sumDemands = (sum(omega_d,1))'
	sumDemandDeviations = sumDemands - sumMeanDemands*ones(Float64,numScenarios)

	# compute the maximum wind power outputs for the scenarios
	pmax_scen = pmax*ones(Float64,1,numScenarios)
	pmax_scen[numRegGen+1:end,:] = omega_w

	# initialize model
	if(useGurobi)
		mod = Model(solver=GurobiSolver(gurobi_env,Presolve=0,OutputFlag=0,BarConvTol=1E-06,OptimalityTol=1E-06))
	else
		mod = Model(solver=IpoptSolver(tol=1E-06,print_level=0))
	end

	# first-stage variables
	@variable(mod, pmin[i] <= p0[i=1:numGen] <= pmax_FS[i])
	@variable(mod, 0 <= rp[i=1:numGen] <= pmax_FS[i] - pmin[i])
	@variable(mod, 0 <= rm[i=1:numGen] <= pmax_FS[i] - pmin[i])
	@variable(mod, 0 <= alpha[i=1:numGen] <= 1)

	# recourse variables
	@variable(mod, p[1:numGen,1:numScenarios])
	@variable(mod, theta[1:numBuses,1:numScenarios])

	# auxiliary variables
	@variable(mod, posResViol[1:numGen,1:numScenarios] >= 0)
	@variable(mod, negResViol[1:numGen,1:numScenarios] >= 0)
	@variable(mod, lineViol[1:numLines,1:numScenarios] >= 0)
	@variable(mod, genBoundViol[1:numGen,1:numScenarios] >= 0)

	# objective function
	@objective(mod, Min, sum(genCost[i]*p0[i] + posResCost[i]*rp[i] + negResCost[i]*rm[i] for i=1:numGen) + (1.0/numScenarios)*sum(gamma_res*posResCost[i]*posResViol[i,scen] + gamma_res*negResCost[i]*negResViol[i,scen] for scen=1:numScenarios for i = 1:numGen) + (1.0/numScenarios)*sum(gamma_line*(lineViol[k,scen])^2 for k = 1:numLines for scen = 1:numScenarios) + (1.0/numScenarios)*sum(gamma_gen*(genBoundViol[i,scen])^2 for i = 1:numGen for scen = 1:numScenarios))

	# first-stage constraints. same as in the reserve saturation model
	@constraint(mod, sum(p0[i] for i=1:numGen) == sumMeanDemands)
	@constraint(mod, [i=1:numRegGen], p0[i] + rp[i] <= pmax[i])
	@constraint(mod, [i=1:numRegGen], p0[i] - rm[i] >= pmin[i])
	@constraint(mod, [i=1:numGen; genGuarRes[i]], alpha[i] >= min_alpha)
	@constraint(mod, [i=1:numRegGen; !genRes[i]], alpha[i] == 0)
	@constraint(mod, sum(alpha[i] for i=1:numGen) == 1)

	# for the no wind reserves case, set the wind generator participation factor to be small enough
	if(caseID == "no_wind_reserves")
		@constraint(mod, [i=1:numWindGen], alpha[i+numRegGen] == min_alpha/10.0)
	end

	# constraints defining p: affine recourse policy without worrying about generation limits
	@constraint(mod, [i=1:numGen,scen=1:numScenarios], p[i,scen] == p0[i] + alpha[i]*sumDemandDeviations[scen])

	# phase angle constraints
	@constraint(mod, [i=1:numBuses,scen=1:numScenarios], sum(betaAdjBus[i][k]*(theta[i,scen] - theta[linesAdjBus[i][k],scen]) for k = 1:size(linesAdjBus[i],1)) == sum(p[regGenAtBus[i][k],scen] for k = 1:size(regGenAtBus[i],1)) + sum(p[windGenAtBus[i][k],scen] for k = 1:size(windGenAtBus[i],1)) - sum(omega_d[loadAtBus[i][k],scen] for k = 1:size(loadAtBus[i],1)))
	# reference bus
	@constraint(mod, [scen=1:numScenarios], theta[numBuses,scen] == 0)

	# line violations
	@constraint(mod, [k=1:numLines,scen=1:numScenarios], lineViol[k,scen] >= beta[k]*(theta[lines[k][1],scen] - theta[lines[k][2],scen]) - delta_line*flow_max[k])
	@constraint(mod, [k=1:numLines,scen=1:numScenarios], lineViol[k,scen] >= -beta[k]*(theta[lines[k][1],scen] - theta[lines[k][2],scen]) - delta_line*flow_max[k])

	# auxiliary constraints for violations of positive and negative reserves
	@constraint(mod, [i=1:numGen,scen=1:numScenarios], posResViol[i,scen] >= p[i,scen] - p0[i] - rp[i])
	@constraint(mod, [i=1:numGen,scen=1:numScenarios], negResViol[i,scen] >= p0[i] - p[i,scen] - rm[i])

	# generator bounds violation
	@constraint(mod, [i=1:numGen,scen=1:numScenarios], genBoundViol[i,scen] >= p[i,scen] - (pmin[i] + delta_gen*(pmax_scen[i,scen]-pmin[i])))
	@constraint(mod, [i=1:numGen,scen=1:numScenarios], genBoundViol[i,scen] >= (pmin[i] + (1-delta_gen)*(pmax_scen[i,scen]-pmin[i])) - p[i,scen])

	solve_status = solve(mod)

	if(solve_status != :Optimal && solve_status != :Suboptimal)
		if(useGurobi)
			println("GP model not solved to (sub)optimality using Gurobi!!!")
		else
			println("GP model not solved to (sub)optimality using Ipopt!!!")
		end
		println("solve_status: ",solve_status)
	end

	p0_soln = getvalue(p0)
	rp_soln = getvalue(rp)
	rm_soln = getvalue(rm)
	alpha_soln = getvalue(alpha)

	return p0_soln, rp_soln, rm_soln, alpha_soln, solve_status
end

#*=======================================
#*============ RUNS START HERE ==========
#*=======================================

tic()

# number of line flow violation penalty iterations
const numGammaLineIter = size(gamma_line_val,1)
# number of generator bounds violation penalty iterations
const numGammaGenIter = size(gamma_gen_val,1)

# print solutions to the screen at each iteration?
const printIterSoln = false

# directory name for storing results
const baseDirName = "C:/Users/rkannan/Desktop/DCOPF-reserve-saturation/experiments/" * model * "/" * "gp_" * caseID
mkpath(baseDirName)

#*========= PRINT OPTIONS ====================
const storeResults = true

const solnFile = "solution.txt"

const p0File = "p0_soln.txt"
const rpFile = "rp_soln.txt"
const rmFile = "rm_soln.txt"
const alphaFile = "alpha_soln.txt"

const costFile = "expected_cost.txt"
const costPowerFile = "expected_cost_power.txt"
const costReservesFile = "expected_cost_reserves.txt"
const costPosResPenFile = "expected_cost_pos_res_pen.txt"
const costNegResPenFile = "expected_cost_neg_res_pen.txt"

const sumLineViolFile = "sumLineViol.txt"
const lineViolProbFile = "lineViolProb.txt"

const powerGenFile = "expected_power_gen.txt"
const windSpillFile = "expected_wind_spill.txt"

const avgLineViolFile = "avgAbsLineFlowViol.txt"
const genBoundViolProbFile = "genBoundViolProb.txt"
const genViolFile = "genBoundViol.txt"

const timeFile = "solutionTime.txt"

#*========= STORE RESULTS ====================
if(storeResults)

	subDirName = baseDirName * "/"
	mkpath(subDirName)

	# write details to text file, including some key details about the test instance
	details_file = subDirName * model * ".txt"
	open(details_file, "w") do f
		write(f,"Model: $model \n")
		write(f,"caseID: $caseID \n")
		write(f,"Number of scenarios (numScenarios): $numScenarios \n\n")

		write(f,"======================\n")
		write(f,"Generator bounds violation penalty parameters \n")
		write(f,"----------------------\n")
		write(f,"Penalty factors (gamma_gen): $gamma_gen_val \n")
		write(f,"Beyond what fraction to penalize? (delta_gen): $delta_gen \n\n")

		write(f,"======================\n")
		write(f,"Line flow violation penalty parameters \n")
		write(f,"----------------------\n")
		write(f,"Penalty factors (gamma_line): $gamma_line_val \n")
		write(f,"Beyond what fraction to penalize? (delta_line): $delta_line \n\n")

		write(f,"======================\n")
		write(f,"Reserves parameters \n")
		write(f,"----------------------\n")
		write(f,"How much costlier are reserves compared to normal generation? (reservesCostFactor): $reservesCostFactor \n")
		write(f,"What is the penalty factor for exceeding reserves? (gamma_res): $gamma_res \n")
		write(f,"How much costlier are wind reserves compared to the cheapest regular generator? (windReservesCostFactor): $windReservesCostFactor \n")
		write(f,"\n\n")

		write(f,"**********************\n")
		write(f,"SOME PARAMETERS OF THE INPUT MODEL \n\n")
		write(f,"**********************\n")

		write(f,"======================\n")
		write(f,"Generator bounds \n")
		write(f,"----------------------\n")
		write(f,"Regular generator upper limits: $regGenMaxPow \n")
		write(f,"Minimum participation factor: $min_alpha \n\n")

		write(f,"======================\n")
		write(f,"Generator cost parameters \n")
		write(f,"----------------------\n")
		write(f,"Regular generator linear cost coefficients: $regGenCost \n")
		write(f,"Regular generator positive reserves cost coefficients: $regGenPosResCost \n")
		write(f,"Regular generator negative reserves cost coefficients: $regGenNegResCost \n")
		write(f,"Wind generator positive reserves cost coefficients: $windGenPosResCost \n")
		write(f,"Wind generator negative reserves cost coefficients: $windGenNegResCost \n\n")

		write(f,"======================\n")
		write(f,"Line flow parameters \n")
		write(f,"----------------------\n")
		write(f,"Line flow factor: $lineFlowFactor \n")
		write(f,"Line flow limits: $flow_max \n\n")

		write(f,"======================\n")
		write(f,"Demand parameters \n")
		write(f,"----------------------\n")
		write(f,"Demand factor: $demandFactor \n")
		write(f,"Average loads/demands: $meanDemand \n")
		write(f,"Demand stdev factor: $demandStdevFactor \n")
		write(f,"Standard deviation of loads/demands: $demandStdev \n")
		write(f,"Demand covariance matrix: $demandCovMat \n\n")

		write(f,"======================\n")
		write(f,"Wind parameters \n")
		write(f,"----------------------\n")
		write(f,"Wind factor: $windFactor \n")
		write(f,"Average wind output: $windPredMean \n")
		write(f,"Wind prediction stdev factor: $windPredStdevFactor \n")
		write(f,"Standard deviation of wind output: $windPredStdev \n")
		write(f,"Do the wind generators provide reserves?: $windGenRes \n")
	end

	# write the line penalty coefficients to a file
	gamma_line_file = subDirName * "gamma_line.txt"
	open(gamma_line_file,"w") do f
		for iter_line = 1:numGammaLineIter
			write(f,"$(gamma_line_val[iter_line])\n")
		end
	end

	# write the generator penalty coefficients to a file
	gamma_gen_file = subDirName * "gamma_gen.txt"
	open(gamma_gen_file,"w") do f
		for iter_gen = 1:numGammaGenIter
			write(f,"$(gamma_gen_val[iter_gen])\n")
		end
	end

end
#*============================================

# replicate the experiments multiple times, if needed
for rep = 1:numReplicates

	println("\n******************************")
	println("******************************")
	println("REPLICATE ",rep,"/",numReplicates)
	println("******************************")

	tic()

	#*============ SCENARIO GENERATION ==========
	# for each replicate, generate common set of scenarios that are used across the various penalty settings
	omega_d, omega_w = generateScenarios(numScenarios)

# in case we want to solve for more replicates later
	if(rep >= startingRep)

		# sweep through the line limit penalty parameters
		for iter_line = 1:numGammaLineIter

			#*========= STORE RESULTS ====================
			if(storeResults)

				subDirName2 = baseDirName * "/rep" * string(rep) * "/lp" * string(iter_line) * "/"
				mkpath(subDirName2)

				cost_file = subDirName2 * costFile
				open(cost_file, "w") do f
				end

				cost_power_file = subDirName2 * costPowerFile
				open(cost_power_file, "w") do f
				end

				cost_reserves_file = subDirName2 * costReservesFile
				open(cost_reserves_file, "w") do f
				end

				cost_pos_res_file = subDirName2 * costPosResPenFile
				open(cost_pos_res_file, "w") do f
				end

				cost_neg_res_file = subDirName2 * costNegResPenFile
				open(cost_neg_res_file, "w") do f
				end

				sum_viol_file = subDirName2 * sumLineViolFile
				open(sum_viol_file, "w") do f
				end

				viol_prob_file = subDirName2 * lineViolProbFile
				open(viol_prob_file, "w") do f
				end

				power_gen_file = subDirName2 * powerGenFile
				open(power_gen_file, "w") do f
				end

				wind_spill_file = subDirName2 * windSpillFile
				open(wind_spill_file, "w") do f
				end

				avg_viol_file = subDirName2 * avgLineViolFile
				open(avg_viol_file, "w") do f
					write(f,"Average/Maximum line flow violations: \n\n")
				end

				gen_viol_prob_file = subDirName2 * genBoundViolProbFile
				open(gen_viol_prob_file, "w") do f
				end

				gen_viol_file = subDirName2 * genViolFile
				open(gen_viol_file, "w") do f
					write(f,"Probability/Median/Maximum generator bound violations: \n\n")
				end

				soln_file = subDirName2 * solnFile
				open(soln_file, "w") do f
				end

				p0_file = subDirName2 * p0File
				open(p0_file, "w") do f
				end

				rp_file = subDirName2 * rpFile
				open(rp_file, "w") do f
				end

				rm_file = subDirName2 * rmFile
				open(rm_file, "w") do f
				end

				alpha_file = subDirName2 * alphaFile
				open(alpha_file, "w") do f
				end

				time_file = subDirName2 * timeFile
				open(time_file, "w") do f
				end

			end
			#*============================================

			println("\n==============================")
			println("==============================")
			println("LINE PENALTY ITER ",iter_line,"/",numGammaLineIter," with penalty ",gamma_line_val[iter_line])
			println("==============================")

			tic()

			# sweep through the penalty coefficients on generator bounds violation
			for iter_gen = 1:numGammaGenIter

				println("\n*********************")
				println("GEN PENALTY ITER ",iter_gen,"/",numGammaGenIter," with penalty ",gamma_gen_val[iter_gen],"\n")

				tic()

				p0_soln, rp_soln, rm_soln, alpha_soln, solve_status = solveGPModel(gamma_gen_val[iter_gen],delta_gen,gamma_line_val[iter_line],delta_line,omega_d,omega_w)

				# if Gurobi did not solve the model as expected, try IPOPT
				if(solve_status != :Optimal && solve_status != :Suboptimal)
					println("--- Gurobi did not converge! Going to try IPOPT ---")
					p0_soln, rp_soln, rm_soln, alpha_soln, solve_status = solveGPModel(gamma_gen_val[iter_gen],delta_gen,gamma_line_val[iter_line],delta_line,omega_d,omega_w,false)
				end

				solnTime::Float64 = toq()

				if(solve_status != :Optimal && solve_status != :Suboptimal)
					println("*****SKIPPING ITERATION BECAUSE BOTH GUROBI AND IPOPT DID NOT CONVERGE!*****")
				else
					p0_soln, rp_soln, rm_soln, alpha_soln, ~ = projectOntoFeasibleSet(p0_soln, rp_soln, rm_soln, alpha_soln)

					if(printIterSoln)
						println("GP solution. Regular generator variables first \n")
						println("gamma_line: ",gamma_line_val[iter_line],"  gamma_gen: ",gamma_gen_val[iter_gen])
						println("p0: ",round.(p0_soln,displayPrecision))
						println("rp: ",round.(rp_soln,displayPrecision))
						println("rm: ",round.(rm_soln,displayPrecision))
						println("alpha: ",round.(alpha_soln,displayPrecision))
						println("")
					end

					tic()

					# estimate quality of solution
					cost_soln::Float64, cost_power_soln::Float64, lineViolProb_soln::Float64, sumLineViol_soln::Float64, cost_reserves_soln::Float64, cost_pos_res_soln::Float64, cost_neg_res_soln::Float64, power_gen_soln, wind_spill_soln, ~, avgAbsLineFlowViol_soln, maxAbsLineFlowViol_soln = estimateCostOfSoln(p0_soln, rp_soln, rm_soln,alpha_soln,tau_sat,gamma_res,tau_pos,gamma_line_val[iter_line],delta_line,true,true)

					probAnyGenBoundViol_soln::Float64, probGenBoundViol_soln, maxGenBoundViol_soln, medianGenBoundViol_soln = estimateGenBoundViolStats(p0_soln,alpha_soln,true,false)

					costTime::Float64 = toq()

					totalTimeIter::Float64 = solnTime + costTime

					println("ITER TIME: ",round.(totalTimeIter,2),"  SOLN TIME: ",round.(solnTime,2),"  COST TIME: ",round.(costTime,2),"\n")

					#*========= STORE RESULTS ==========
					if(storeResults)

						open(cost_file, "a") do f
							write(f,"$cost_soln \n")
						end

						open(cost_power_file, "a") do f
							write(f,"$cost_power_soln \n")
						end

						open(cost_reserves_file, "a") do f
							write(f,"$cost_reserves_soln \n")
						end

						open(cost_pos_res_file, "a") do f
							write(f,"$cost_pos_res_soln \n")
						end

						open(cost_neg_res_file, "a") do f
							write(f,"$cost_neg_res_soln \n")
						end

						open(sum_viol_file, "a") do f
							write(f,"$sumLineViol_soln \n")
						end

						open(viol_prob_file, "a") do f
							write(f,"$lineViolProb_soln \n")
						end

						open(power_gen_file, "a") do f
							for i = 1:numGen
								write(f,"$(power_gen_soln[i]) ")
							end
							write(f,"\n")
						end

						open(wind_spill_file, "a") do f
							for i = 1:numWindGen
								write(f,"$(wind_spill_soln[i]) ")
							end
							write(f,"\n")
						end

						open(avg_viol_file, "a") do f
							write(f,"ITER NUM. $iter_gen \n")
							for k = 1:numLines
								write(f,"($(lines[k][1]),$(lines[k][2])): $(avgAbsLineFlowViol_soln[k])/$(maxAbsLineFlowViol_soln[k])  \n")
							end
							write(f,"\n")
						end

						open(gen_viol_prob_file, "a") do f
							write(f,"$probAnyGenBoundViol_soln \n")
						end

						open(gen_viol_file, "a") do f
							write(f,"ITER NUM. $iter_gen \n")
							for i = 1:numGen
								write(f,"Gen #$i: $(probGenBoundViol_soln[i])/$(medianGenBoundViol_soln[i])/$(maxGenBoundViol_soln[i]) \n")
							end
							write(f,"\n")
						end

						open(soln_file, "a") do f
							write(f,"ITER NUM. $iter_gen \n")
							write(f,"Expected cost: $cost_soln \n")
							write(f,"Expected cost without line penalty costs: $cost_power_soln \n")
							write(f,"Expected penalty for exceeding positive reserves: $cost_pos_res_soln \n")
							write(f,"Expected penalty for exceeding negative reserves: $cost_neg_res_soln \n")
							write(f,"Expected sum of line flow limit violations: $sumLineViol_soln \n")
							write(f,"ANY line limit violation prob: $lineViolProb_soln \n")
							write(f,"ANY generator bounds violation prob: $probAnyGenBoundViol_soln \n")
							write(f,"p0: $p0_soln \n")
							write(f,"rp: $rp_soln \n")
							write(f,"rm: $rm_soln \n")
							write(f,"alpha: $alpha_soln \n\n")
							write(f,"Expected power generation: $power_gen_soln \n")
							write(f,"Expected wind spill: $wind_spill_soln \n")
							write(f,"\n\n")
						end

						open(p0_file, "a") do f
							for i = 1:numGen
								write(f,"$(p0_soln[i]) ")
							end
							write(f,"\n")
						end

						open(rp_file, "a") do f
							for i = 1:numGen
								write(f,"$(rp_soln[i]) ")
							end
							write(f,"\n")
						end

						open(rm_file, "a") do f
							for i = 1:numGen
								write(f,"$(rm_soln[i]) ")
							end
							write(f,"\n")
						end

						open(alpha_file, "a") do f
							for i = 1:numGen
								write(f,"$(alpha_soln[i]) ")
							end
							write(f,"\n")
						end

						open(time_file, "a") do f
							write(f,"$totalTimeIter \n")
						end
					end
					#*==================================

				end # solve_status

			end # generator penalty iterations

			totalTime_1::Float64 = toq()

			println("Total time for line penalty iteration #",iter_line,": ",round.(totalTime_1,2))

			#*========= STORE RESULTS ==========
			if(storeResults)
				open(time_file, "a") do f
					write(f,"\n\nTOTAL TIME: $totalTime_1 \n")
				end
			end
			#*==================================

		end # line penalty iterations

	end # end replicate num check

	replicateTime::Float64 = toq()
	println("Time for replicate #",rep,":  ",round.(replicateTime,2))

end # replicates

toc() # total time for all runs

max_stars_repo_name: rohitkannan/DCOPF-reserve-saturation

max_issues_repo_name: rohitkannan/DCOPF-reserve-saturation

max_forks_repo_name: rohitkannan/DCOPF-reserve-saturation


text:
function distance(s1::AbstractString, s2::AbstractString)
    length(s1) != length(s2) && throw(ArgumentError("Sequences must have the same length"))
    reduce(+, a != b for (a, b) in zip(s1, s2); init=0)
end


text:
using UUIDs

"""
`AbstractJobStore` is an abstract type for jobstores

A jobstore is a data structure which is responsible of
storing jobs that should be executed later.
"""
abstract type AbstractJobStore end

"""
    get_job_id()

Returns a job identifier (a `job_id`).

It's preferable to use `get_job_id(jobstore)` to ensure
that a `job_id` is unique for a given `JobStore`.
"""
function get_job_id()
    string(UUIDs.uuid4())
end

include("memory.jl")


text:
# This files implements Lovasz Theta computation described in the following paper:
#
# An SDP Primal-Dual Algorithm for Approximating the Lovász-Theta Function
# T.-H. Hubert Chan, Kevin L. Chang and Rajiv Raman
# Algorithmica (2014) 69:605-618
# DOI: 10.1007/s00453-013-9756-5
#
using Base.Test
Pkg.add("Graphs")
Pkg.add("SVM")

using SVM

function test_expm()
  a = 2e10
  b = 4e8/6
  c = 200/3
  d = 3
  e = 1e-8
  A = [0 e 0; -(a+b) -d a; c 0 -c]
  E = expm(A)
  E_true = [  0.446849 1.54044e-9 0.462811; -5.74307e6 -0.015283 -4.52654e6;    0.447723 1.5427e-9 0.463481]
  Z = E_true ./ E # test relative error
  @test_approx_eq_eps ones(3,3) Z 1e-3
end
test_expm()

function ORACLE(W, β, E)
  n = size(W, 1)
  Ø = zeros(n, n)
  if isequal(W, Ø)
    z = β
    Y = Ø
    return ("dual", z, Y)
  else
    J = ones(n, n)
    t = trace(W)
    X = W ./ t
    σ = sum (J .* W) - β
    if σ <= 0.0
      z = β
      Y = Ø
      return ("dual", z, Y)
    end
    X_E = X .* E
    f = normfro(X_E)
    if f >= σ/(β*n)
      z = β
      Y = σ/(f^2) * X_E
      return ("dual", z, Y)
    end
    X_tilde = X - X_E + σ/(β*n)*J
    X_hat  = X_tilde ./ trace(X_tilde)
    return ("primal", null, X_hat)
  end
end

function ϑ(E, ρ, δ, β)
  n = size(E, 1)
  J = ones(n, n)
  I = eye(n)
  S_M = zeros(n, n)
  S_Y = zeros(n, n)
  W = I
  ϵ = δ/(2*ρ)
  ϵ_1 = -log(1 - ϵ)
  T = iceil ( 8 * ρ^2 * log(n) /  δ^2 )
  for t=1:T
    (solution, val, mat) = ORACLE(W, β, E)
    if solution == "primal"
      return (solution, val, mat)
    end
    Y = sum( mat .* E)
    M = ( β*I + Y - J + ρ*I ) / (2*ρ)
    S_M += M
    S_Y += Y
    print(W)
    W = expm(-ϵ_1 * S_M)
  end
  z = β + δ
  Y_bar = 1\T * S_Y
  return ("dual", z, Y_bar)
end

using Graphs

# function test_ϑ()
  n = 10
  p = 0.5
  G = Graphs.erdos_renyi_graph(n, p)
  E = int(Graphs.adjacency_matrix(G))
  β = n
  ρ = β*n
  δ = 0.5
  W = eye(n, n)
  (sol, val, mat) = ORACLE(W, β, E)
  print(sol)
(solution, val, mat) = ϑ(E, ρ, δ, β)
# end

test_ϑ()