text:
"""
	writeatmosph(data,fileout;decimal)
Write DataFrame to Hydrus1D Atmospheric data format

**Input**
* data: dataframe containing columns + timevector in DateTime format.
* fileout: full name of the output file
* decimal: optional output precision for either all or individual columns.
> input dataframe should contain :Prec,:rSoil,:rRoot,:hCritA,:rB,:hB,:ht columns,
> if not, will be set to zero. Tested for  Hydrus-1D 4.16.0110

**Example**
```
dataout = DataFrame(Prec=[0.01,0.1,0.2,0.3],
				 rSoil=[0.0,0.1,0.2,0.9],
	   			 hCritA=[1,1,1,1]
	   			 datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,1),
		   			  DateTime(2010,1,1,2),DateTime(2010,1,1,3)],);
writeatmosph(dataout,"test/output/atmosph_data.in",
			decimal=[1],
			hCritS=1000.,Extinction=-1.0);
```
"""
function writeatmosph(data::DataFrame,fileout::String;decimal=[4],hCritS=0.0,Extinction=-1.0)
	channels,timei = FileTools.findchannels(data);
	# Round data to output precision
	dataout = FileTools.round2output(data,decimal,channels);
	open(fileout,"w") do fid
		# Write header
		writeatmosph_head(fid,size(data,1),Extinction,hCritS);
		# Write data
		for i = 1:size(dataout,1)
			# Write date-index
			@printf(fid,"%11g",i);
			# add data
			for j in [:Prec,:rSoil,:rRoot,:hCritA,:rB,:hB,:ht]
				if haskey(dataout,j)
					@printf(fid,"%12g",dataout[j][i]);
				else # set missing columns to zero
					@printf(fid,"%12g",0);
				end
			end
			@printf(fid,"\n");
		end
		# add footer
		@printf(fid,"end*** END OF INPUT FILE 'ATMOSPH.IN' **********************************\n");
	end
end

"""
auxiliary function to write fixed header
"""
function writeatmosph_head(fid::IOStream,l::Int,extinc::Float64,h::Float64)
	@printf(fid,"Pcp_File_Version=4\n");
	@printf(fid,"*** BLOCK I: ATMOSPHERIC INFORMATION  ");
	for i in 1:34 @printf(fid,"*");end
	@printf(fid,"\n   MaxAL                    (MaxAL = number of atmospheric data-records)\n")
	@printf(fid,"   %i\n",l);
	@printf(fid," DailyVar  SinusVar  lLay  lBCCycles lInterc lDummy  lDummy  lDummy  lDummy  lDummy\n");
	if extinc > -1e-6
		@printf(fid,"       f       f       t       f       f       f       f       f       f       f\n")
		@printf(fid," Extinction\n   %.6g\n",extinc);
	else
		@printf(fid,"       f       f       f       f       f       f       f       f       f       f\n")
	end
	@printf(fid," hCritS                 (max. allowed pressure head at the soil surface)\n");
	@printf(fid,"   %.6g\n",h);
	@printf(fid,"       tAtm        Prec       rSoil       rRoot      hCritA");
	@printf(fid,"          rB          hB          ht      RootDepth\n");
end

"""
	writeprofile1d(nodeinfo,fileout)
Write Hydrus1D profile (nodal) file. Suitable only for pure water flow simulation

**Input**
soilinfo: dataframe containg information on:
* :start = soil layer starting depth
* :stop = soil layer end depth
* :res = vertical resolution (constant sampling)
* :h = Initial value of the pressure head
* :Mat = index of the material
* :Lay = subregion number assigned to nodes within the soil layer

fileout: output file name
iObs: list (vector) of the observation nodes for which values of the pressure head & the water content are printed

> No information about temperature or concentrations is used! Hydrus-1D 4.16.0110

**Example**
```
soilinfo = DataFrame(start=[0],stop=[10],res=[0.01],h=[100],Mat=[5],Lay=[1])
output_file = "test/output/profile1d_data.in"
print_nodes = [1,2,3];
writeprofile1d(soilinfo,output_file,iObs=print_nodes);
```
"""
function writeprofile1d(soilinfo::DataFrame,fileout::String;iObs=[0])
	open(fileout,"w") do fid
		# write header
		@printf(fid,"    2\n");# fixed number for "fixed" nodes
		@printf(fid,"%5g%15g%15g%15g\n",1,soilinfo[:start][1],1,1);
		@printf(fid,"%5g%15g%15g%15g\n",2,-soilinfo[:stop][end],1,1);
		@printf(fid,"%5g%5g%5g%5g",getnumofnodes(soilinfo),0,0,1);
		@printf(fid," x         h       Mat  Lay      Beta           Axz            ");
		@printf(fid,"Bxz            Dxz            Temp          Conc           SConc\n");
		# Write nodes
		c = 1;
		for i in 1:size(soilinfo,1)
			n = soilinfo[:start][i]:soilinfo[:res][i]:soilinfo[:stop][i];
			for j in n
				@printf(fid,"%5g%15g%15g%5g%5g",
						c,-j,soilinfo[:h][i],soilinfo[:Mat][i],soilinfo[:Lay][i]);
				# add constant values
				@printf(fid,"  0.000000e+000  1.000000e+000  1.000000e+000  1.000000e+000              \n");
				c += 1; # count nodes
			end
		end
		# add footer
		if iObs[1] != 0
			@printf(fid,"%5g\n",length(iObs));
			for i in iObs
				@printf(fid,"%5g",i)
			end
			@printf(fid,"\n");
		else
			@printf(fid,"    0\n");
		end
	end
end

"""
auxiliary function to get total number of observation points
"""
function getnumofnodes(soilinfo)
	c = 0;
	for i in 1:size(soilinfo,1)
		for j in soilinfo[:start][i]:soilinfo[:res][i]:soilinfo[:stop][i]
			c += 1; # count nodes
		end
	end
	return c
end

"""
	obnode = readhydrus1d_obsnode(filein,paramout)
Read Hydrus1D observation node output

**Input**
* filein: input file name
* paramout: what parameter should be returned (theta,h or Flux)

**Output**
obnode: dataframe containg:
* :time = time index
* :paramoutX = parameter value for node number X

> Designed only for Hydrus1D with water fluxes only! Hydrus-1D 4.16.0110

**Example**
```
input_file = "test/input/hydrus1d_Obs_Node.out"
obnode = readhydrus1d_obsnode(input_file,paramout=:theta)
```
"""
function readhydrus1d_obsnode(filein::String;paramout::Symbol=:theta)
	# read all data
	dataall = readdlm(filein,skipstart=11,comments=true,comment_char='e');
	obsnodes = getnodenumber(filein,paramout);
	dataout = DataFrame(time=dataall[:,1]);
	if paramout == :h
		col = 2;
	elseif paramout == :Flux
		col = 4;
	else
		col = 3;
	end
	for (i,v) in enumerate(obsnodes)
		dataout[v] = dataall[:,col+(i-1)*3];
	end
	return dataout;
end

"""
	atm = readatmosph(filein)
Read Hydrus1D ATMOSPH.IN file

**Input**
* filein: input file name

**Output**
obnode: dataframe containg:
* :time = time index
* :param = parameter name

> Designed only for Hydrus1D Hydrus-1D 4.16.0110

**Example**
```
input_file = "test/input/hydrus1d_atmosph.in"
atm = readatmosph(input_file)
```
"""
function readatmosph(filein::String)
	nh = headlines(filein,"tAtm");
	temp = readdlm(filein,skipstart=nh-1);
	out = DataFrame(time=temp[2:end-1,1]);
	for (i,v) in enumerate(temp[1,2:end])
		out[Symbol(v)] = eltype(temp[2,i+1])==Char ?
						zeros(Float64,length(temp[2:end-1,1])).*NaN :
						convert(Vector{Float64},temp[2:end-1,i+1])
	end
	return out
end

"""
Aux function to count number of header lines (up to a given string)

**Example**
```
filein = "f:/mikolaj/code/data_processing/sites/aggo/hydro/hydrus/AGGO/ATMOSPH.IN";
nh = headlines(filein,"tAtm")
```
"""
function headlines(filein::String,flag::String)
	c = 0;
	open(filein,"r") do fid
		row = readline(fid);
		while !occursin(flag,row)
			row = readline(fid);
			c += 1;
			if eof(fid)
				c = 0; break
			end
		end
	end
	return c+1
end

"""
auxiliary function to read node numbers
"""
function getnodenumber(filein::String,paramout::Symbol)
	obsnodes = Vector{Symbol}();
	# read node numbers
	open(filein,"r") do fid
		row = "";
		for i in 1:9
			row = readline(fid);
		end
		obsnodes = parsenodenumber(row,paramout);
	end
	return obsnodes;
end

"""
auxiliary function to extract numbers from string
"""
function parsenodenumber(row::String,paramout::Symbol)
	temp = split(row,"Node(");
	c = 1000; # max obs node number in the hydrusFile is 999 (above just ***)
	out = Vector{Symbol}();
	for i in temp[2:end] # starts with empty field (see split)
		if !occursin("*",i) # numbers below 1000
			obsnode = Meta.parse(split(i,")")[1]) # ends with ')'
		else
			obsnode = c;
			c += 1;
		end
		out = vcat(out,Symbol(string(paramout)*string(obsnode)));
	end
	return out
end

"""
	obnode = readhydrus1d_nedinf(filein,paramout)
Read Hydrus1D Node Info file (Nod_inf.out)

**Input**
* filein: input file name
* paramout: what parameter should be returned (all (default), theta, h or k)

**Output**
theta,k,h: dataframe containg:
* :time = time index
* :nodeX = valus for node number X

> Designed only for Hydrus1D with water fluxes only!

**Example**
```
input_file = "test/input/hydrus1d_Nod_Inf.out"
moisture = readhydrus1d_nedinf(input_file,paramout=:theta)
# moisture,k,h = readhydrus1d_nedinf(input_file)
```
"""
function readhydrus1d_nodinf(filein;paramout::Symbol=:all)
	# declare output
	timeout = -1;
	h = Matrix{Float64}(undef,0,0);
	theta = Matrix{Float64}(undef,0,0);
	k = Matrix{Float64}(undef,0,0);
	node = Vector{Int}();
	depth = Vector{Float64}();
	# read header
	open(filein,"r") do fid
		row = readhydrus1d_nodinf_head(fid);
		while !eof(fid)
			timeblock,datablock = readhydrus1d_nodinf_block(fid,row);
			if typeof(timeout)==Int # initial run
				timeout = timeblock;
				h = datablock[:,3]; # get Heat values (fixed position)
				theta = datablock[:,4];
				k = datablock[:,5]; # get hydraulic conductivity values (fixed position)
				node = round.(Int,datablock[:,1]); # read nodes and depth only once
				depth = datablock[:,2];
			else
				timeout = vcat(timeout,timeblock);
				h = hcat(h,datablock[:,3]);
				theta = hcat(theta,datablock[:,4]);
				k = hcat(k,datablock[:,5]);
			end
			row = readline(fid); # to reach EOF after last node info block
		end
	end
	if paramout == :theta # could not pass "paramout" as input
		return readhydrus1d_nodinf_df(timeout,node,theta)
	elseif paramout == :k
 		return readhydrus1d_nodinf_df(timeout,node,k)
	elseif paramout == :h
 		return readhydrus1d_nodinf_df(timeout,node,h)
	else # return 3 dataframes (default)
		return readhydrus1d_nodinf_df(timeout,node,theta),
				readhydrus1d_nodinf_df(timeout,node,k),
				readhydrus1d_nodinf_df(timeout,node,h)
    end
end

"""
auxiliary function to read time and data for each block
"""
function readhydrus1d_nodinf_block(fid,row)
	# find time stamp
	while !occursin("Time:",row)
		row = readline(fid);
	end
	timeout = Meta.parse(split(row,":")[2])
	# ignore header for each time step
	row = readhydrus1d_nodinf_head(fid);
	# read data
	datanode = readhydrus1d_nodinf_data(fid);
	return timeout,datanode
end

"""
auxiliary function to read file header
"""
function readhydrus1d_nodinf_head(fid)
	row = " ";
	for i in 1:5
		row = readline(fid);
	end
	return row;
end

"""
auxiliary function to get data for current block
"""
function readhydrus1d_nodinf_data(fid)
	# stack the data to a matrix
	datanode = Matrix{Float64}(undef,1,11);
	row = readline(fid);
	while !occursin("end",row)
		datanode = vcat(datanode,row |> split |> x -> Base.parse.(Float64,x) |> transpose);
		row = readline(fid);
	end
	return datanode[2:end,:]; # first line is a dummy used for declaration
end

"""
auxiliary function to convert matrix to dataframe
"""
function readhydrus1d_nodinf_df(timeout,nodes,dataout)
	df = DataFrame(time = timeout);
	for (i,v) in enumerate(nodes)
		temp = Symbol("node"*string(v));
		df[temp] = dataout[i,:];
	end
	return df
end

"""
	corr_hydrus1d_obsnode(file_obs_node,file_profile)
Correct Hydrus1D Observation node file if this does not contain all node numbers
in the header. This is motly the case if more then 10 observation nodes are
selected (manually) in the profile settings.

**Input**
* file_obs_node: Observation node file (mostly "Obs_Node.out"). The header of this file will be overwritten
* file_profile: profile settings file (mostly "PROFILE.DAT"). Will be used to get full row of Node IDs

**Example**
```
file_obs_node = "f:/mikolaj/code/data_processing/sites/aggo/hydro/model/hydrus/AGGO/Obs_Node.out"
file_profile = "f:/mikolaj/code/data_processing/sites/aggo/hydro/model/hydrus/AGGO/PROFILE.DAT"
FileTools.corr_hydrus1d_obsnode(file_obs_node,file_profile)
```
"""
function corr_hydrus1d_obsnode(file_obs_node::String,file_profile::String)
	fileout = "temp_Obs_Node.out";

	## Get node numbers
	temp = ""
	fid = open(file_profile,"r");
	node_id = fid |> readlines |> x->x[end] |> split |> x-> Base.parse.(Int,x)
	close(fid)

	## Write with new Node Info
	open(file_obs_node,"r") do fid
		fid_out = open(fileout,"w");
		while !eof(fid)
			row = readline(fid);
			if occursin("Node(",row)
				for i in node_id
					@printf(fid_out,"                      Node(%3i)",i)
				end
			else
				@printf(fid_out,"%s",row);
			end
			@printf(fid_out,"\n");
		end
		close(fid_out)
	end
	cp(fileout,file_obs_node,force=true);
	rm(fileout);
end


text:
using BellScenario

function profile_test()
    vertices = map( v -> convert.(Int64, v), filter(v -> !in(v, LocalPolytope.vertices((5,1),(1,5))), LocalPolytope.vertices((5,1),(1,5), dits=2)))
    scenario = LocalSignaling(5,5,2)
    println(length(vertices))

    BG = BellGame([1 0 0 0 0;0 1 0 0 0;0 0 1 0 0;0 0 0 1 0;0 0 0 0 1], 2)

    # skip = [BellGame([1 0 0 0 0;1 0 0 0 0;1 0 0 0 0;1 0 0 0 0;0 0 0 0 0], 1)]
    @time facets = LocalPolytope.adjacency_decomposition(vertices, BG, scenario)#,  skip_games=skip)
    println(length(facets))
end

using ProfileView
profile_test()  # run once to trigger compilation (ignore this one)
@profview profile_test()

ext: jl


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIB = MOI.Bridges
const MOIT = MOI.Test

const BRIDGED = MOIB.full_bridge_optimizer(SCIP.Optimizer(display_verblevel=0), Float64)
const CONFIG_BRIDGED = MOIT.Config(atol=5e-3, rtol=1e-4, exclude=Any[
    MOI.ConstraintDual, MOI.ConstraintName, MOI.DualObjectiveValue, MOI.VariableBasisStatus, MOI.ConstraintBasisStatus,
])

@testset "MOI unit tests" begin
    excluded = copy(MOI_BASE_EXCLUDED)
    append!(excluded, [
        "test_linear_Interval_inactive",
        # Can not delete variable while model contains constraints
        "test_linear_integration",
        "test_basic_ScalarQuadraticFunction_ZeroOne",
        "test_basic_VectorAffineFunction_NormOneCone",
        "test_basic_VectorAffineFunction_SOS1",
        "test_basic_VectorAffineFunction_SOS2",
        "test_basic_VectorQuadraticFunction_NormOneCone",
        "test_basic_VectorQuadraticFunction_SOS1",
        "test_basic_VectorQuadraticFunction_SOS2",
        "test_model_ordered_indices", # TODO should fix? ListOf in order of creation
        "test_quadratic_duplicate_terms", # Can not delete variable while model contains constraints
        "test_quadratic_nonhomogeneous", # unsupported by bridge
        "ScalarAffineFunction_ZeroOne",
        "test_model_delete",
        "test_conic_GeometricMeanCone_Vector",
        "test_basic_VectorOfVariables_GeometricMeanCone",
        "test_basic_VectorQuadraticFunction_GeometricMeanCone",
        "test_basic_VectorAffineFunction_GeometricMeanCone",
        "test_variable_delete_SecondOrderCone",
        "test_modification_func_scalaraffine_",
        "test_modification_func_vectoraffine_",
    ])
    MOIT.runtests(
        BRIDGED,
        CONFIG_BRIDGED,
        warn_unsupported=false,
        exclude = excluded,
    )
end


text:
type ClusterInfo
	meta::Dict
	locked::Bool
	cloudname::AbstractString
	skipticks::Bool
	version::AbstractString
	cloudsize::Int
	healthy::Bool
	badnodes::Int
	excludefields::AbstractString
	nodes::Vector
	clouduptimemillis::Int
	nodeidx::Int
	consensus::Bool
	isclient::Bool
end


text:
using Flight.Geodesy
using Flight.Kinematics
using Flight.Atmosphere
using Flight.Airdata
using Flight.Utils

function test_airdata()

    atm_sys = System(AtmosphereDescriptor())
    p = Geographic(alt = AltOrth(13000))
    x_kin = init_x0(KinLTF(), KinInit(v_eOb_b = [300, 10, 30], Ob = p, q_nb = Ry(0.01)))
    dx_pos = copy(x_kin.pos)
    kin_data = f_kin!(dx_pos, x_kin)

    atm_sys.u.wind.v_ew_n[1] = 50
    atm_sys.u.isa_.T_sl += 10 #ISA+10 day

    b = @benchmarkable AirData($kin_data, $atm_sys); run(b)

end


text:
using Documenter, DocumenterTools, COSMO, Literate

@info "Building example problems..."

# utility function from https://github.com/JuliaOpt/Convex.jl/blob/master/docs/make.jl
fix_math_md(content) = replace(content, r"\$\$(.*?)\$\$"s => s"```math\1```")
fix_suffix(filename) = replace(filename, ".jl" => ".md")
function postprocess(cont)
      """
      The source files for all examples can be found in [/examples](https://github.com/oxfordcontrol/COSMO.jl/tree/master/examples/).
      """ * cont
end

# find all example source files
exclude_files = ["chordal_decomposition.jl"; "chordal_decomposition_generate_data.jl"; "maxEigenvalue.jl"; "sum_of_squares.jl"; "lovasz.jl"; "svm_dual.jl"; "svm_dual_simple"; "max_k_cut.jl"];
example_path = joinpath(@__DIR__, "../examples/")
build_path =  joinpath(@__DIR__, "src", "examples/")
files = readdir(example_path)
filter!(x -> endswith(x, ".jl"), files)
filter!(x -> !in(x, exclude_files), files)

for file in files
      Literate.markdown(example_path * file, build_path; preprocess = fix_math_md, postprocess = postprocess, documenter = true, credit = true)
end

# copy some .csv file for the logistic regression example to the build directory
cp(joinpath(@__DIR__, "../examples/chip_data.txt"),
   joinpath(@__DIR__, "src/examples/chip_data.txt"); force = true)

examples_nav = fix_suffix.("./examples/" .* files)

# find all other documentation source files that are build with Literate
example_path = joinpath(@__DIR__, "src", "literate/")
build_path =  joinpath(@__DIR__, "src", "literate", "build/")
files = readdir(example_path)
filter!(x -> endswith(x, ".jl"), files)
for file in files
      Literate.markdown(example_path * file, build_path; preprocess = fix_math_md, documenter = true, credit = true)
end

@info "Makeing documentation..."
makedocs(
  sitename="COSMO.jl",
  authors = "Michael Garstka and contributors.",
  format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        canonical = "https://oxfordcontrol.github.io/COSMO.jl/stable/",
        assets = ["assets/favicon.ico"; "assets/github_buttons.js"; "assets/custom.css"],
        analytics = "UA-134239283-1",
  ),
  pages = [
        "Home" => "index.md",
        "User Guide" => Any[
        "Getting Started" =>  "getting_started.md",
        "JuMP Interface" => "jump.md",
        "Linear System Solver" => "lin_solver.md",
        "Acceleration" => "acceleration.md",
        "Custom Cone Constraint" => "./literate/build/custom_cone.md",
        "Chordal Decomposition" => "decomposition.md",
        "Model Updates" => "./literate/build/portfolio_model_updates.md",
        "Arbitrary Precision" => "./literate/build/arbitrary_precision.md",
        "Performance Tips" => "performance.md"
        ],
        "Method" => "method.md",
        "Examples" => examples_nav,
        "Citing COSMO" => "citing.md",
        "Contributing" => "contributing.md",
        "API Reference" => "api.md"
    ]
)

deploydocs(
    repo = "github.com/oxfordcontrol/COSMO.jl.git")


text:
# fractals.jl

###############################################################################
## The Mandelbrot set
###############################################################################
export Mandelbrot,mandelbrot
# TODO: fractals should be subsets of the complex plane
struct Mandelbrot{T} <: EuclideanDomain{2,T}
    maxiter     ::  Int
    threshold   ::  T
    maskcache   ::  Dict
    box         ::  Domain
end

function mandelbrot(maxiter = 1000, threshold = 1000.0)
    box = rectangle(-1.0, 0.35, -0.65, 0.65)
    M1 = 136
    M2 = 200
    mask1 = computemandelbrotgrid(EquispacedGrid(M1,-1.0, 0.3)× EquispacedGrid(M2,-0.65, 0.65), maxiter, threshold)
    mask2 = computemandelbrotgrid(EquispacedGrid(M1,-1.0, 0.3)× EquispacedGrid(M2,-0.65, 0.65), maxiter, threshold)
    cache = Dict{Int,BitArray{2}}()
    cache[M1] = mask1
    cache[M2] = mask2
    Mandelbrot(maxiter, threshold, cache, box)
end

function mandelbrotiteration(x, maxiter, threshold)
    c = 2*(x[1]+1im*x[2])
    T = typeof(c)
    z = zero(T)
    iter = 0
    while abs(z) < threshold && iter < maxiter
        z = z^2 + c
        iter += 1
    end
    abs(z) < threshold
end

function computemandelbrotgrid(grid, maxiter, threshold)
    mask = BitArray(undef, size(grid))
    fill!(mask,0)
    for i_2 = 1:size(grid, 2)
        for i_1 = 1:size(grid, 1)
            mask[i_1,i_2] = mandelbrotiteration(grid[i_1,i_2], maxiter, threshold)
        end
    end
    mask
end

function indomain_broadcast(grid::ProductGrid, m::Mandelbrot)
    if (infimum(support(grid)) ≈ minimum(boundingbox(m))) && (supremum(support(grid)) ≈ maximum(boundingbox(m)))
        if haskey(m.maskcache, size(grid,1))
            mask = m.maskcache[size(grid,1)]
        else # compute mask and cache it
            mask = computemandelbrotgrid(grid, m.maxiter, m.threshold)
            m.maskcache[size(grid,1)] = mask
        end
    else # Don't cache if the grid doesn't match the bounding box
        mask = computemandelbrotgrid(grid, m.maxiter, m.threshold)
    end
    mask
end

function isapprox(t::Tuple{T,T}, v::SVector{2,T}) where {T}
    return t[1]≈v[1] && t[2]≈v[2]
end
function isapprox(v::SVector{2,T}, t::Tuple{T,T}) where {T}
    return t[1]≈v[1] && t[2]≈v[2]
end
boundingbox(m::Mandelbrot) = m.box

indomain(x::SVector{2}, m::Mandelbrot) = mandelbrotiteration(x, m.maxiter, m.threshold)

show(io::IO, m::Mandelbrot) = print(io, "The Mandelbrot set")

################
## Julia sets
################
export JuliaSet, juliaset
struct JuliaSet{T} <: EuclideanDomain{2,T}
    c           ::  Complex{T}
    maxiter     ::  Int
    maskcache   ::  Dict
    box         ::  Domain
end

function juliaset(c = -0.122565+0.744866im, maxiter = 1000)
    box = rectangle(-0.2, 1.2, -0.4, 0.4)

    mask1 = computejuliasetgrid(EquispacedGrid(200,-0.2, 1.2)× EquispacedGrid(200,-0.4, 0.4), c, maxiter)
    mask2 = computejuliasetgrid(EquispacedGrid(200,-0.2, 1.2)× EquispacedGrid(200,-0.4, 0.4), c, maxiter)
    cache = Dict{Int,BitArray{2}}()
    cache[100] = mask1
    cache[200] = mask2
    JuliaSet(c, maxiter, cache, box)
end

function juliasetiteration(x, c, maxiter)
    gamma = 1 + sqrt(1-4*c)
    z = x[1] + 1im*x[2]
    for i = 1:maxiter
        z = gamma*z*(1-z)
    end
    abs(z) < 1000
end

function computejuliasetgrid(grid, c, maxiter)
    m = size(grid)
    mask = BitArray(undef, m)
    fill!(mask, 0)
    for i_2 = 1:m[2]
        for i_1 = 1:m[1]
            a,b = grid[i_1,i_2]
            mask[i_1,i_2] = juliasetiteration(grid[i_1,i_2], c, maxiter)
        end
    end
    mask
end

function indomain_broadcast(grid, js::JuliaSet)
    if (map(x->infimum(support(x)), elements(grid)) ≈ infimum(boundingbox(js))) && (map(x->supremum(support(x)),elements(grid)) ≈ supremum(boundingbox(js)))
        if haskey(js.maskcache, size(grid,1))
            mask = js.maskcache[size(grid,1)]
        else # compute mask and cache it
            mask = computejuliasetgrid(grid, js.c, js.maxiter)
            js.maskcache[size(grid,1)] = mask
        end
    else # Don't cache if the grid doesn't match the bounding box
        mask = computejuliasetgrid(grid, js.c, js.maxiter)
    end
    mask
end

indomain(x::SVector{2}, js::JuliaSet) = juliasetiteration(x, js.c, js.maxiter)

show(io::IO, js::JuliaSet) = print(io, "A particular Julia Set also known as the Douady rabbit")

boundingbox(js::JuliaSet) = js.box


text:
__precompile__()

module ImplicitSkinning

using SymPy

function grad(f, x)
    [diff(f, x[i]) for i in 1:length(x)]
end

function hrbf_field(x, phi, v, lambda, beta)
    phi_expr = phi(norm(x - v))
    return lambda * phi_expr + dot(beta, grad(phi_expr, x))[1]
end

type HRBFGenerator
	dimension::Integer
	partials::Dict{Tuple{Symbol, Symbol}, Function}
	field_func::Function
end

function HRBFGenerator(dimension::Integer, phi::Function)
	lambda = symbols("lambda", real=true)
	x = Sym[symbols("x$(i)", real=true) for i in 1:dimension]
	v = Sym[symbols("v$(i)", real=true) for i in 1:dimension]
	beta = Sym[symbols("beta$(i)", real=true) for i in 1:dimension]
	f = hrbf_field(x, phi, v, lambda, beta)
	g = grad(f, x)

	df_dlambda_func = lambdify(diff(f, lambda), [x...,v...])
	df_dbeta_func = [lambdify(diff(f, beta[i]), [x...,v...]) for i in 1:dimension]
	dg_dlambda_func = [lambdify(diff(g[i], lambda), [x...,v...]) for i in 1:dimension]
	dg_dbeta_func = [lambdify(diff(g[i], beta[j]), [x...,v...]) for i in 1:dimension, j in 1:dimension]
	f_func = lambdify(f, [lambda, beta..., x..., v...])

	partials = Dict{Tuple{Symbol, Symbol}, Function}()
	partials[(:f, :lambda)] = (point_index, lambda_index) -> (points -> df_dlambda_func(points[point_index,:]..., points[lambda_index,:]...))
	partials[(:f, :beta)] = (point_index, beta_index, coord_index) -> (points -> df_dbeta_func[coord_index](points[point_index,:]..., points[beta_index,:]...))
	partials[(:g, :lambda)] = (point_index, lambda_index, grad_coord_index) -> (points -> dg_dlambda_func[grad_coord_index](points[point_index,:]..., points[lambda_index,:]...))
	partials[(:g, :beta)] = (point_index, beta_index, grad_coord_index, beta_coord_index) -> (points -> dg_dbeta_func[grad_coord_index, beta_coord_index](points[point_index,:]..., points[beta_index,:]...))
	HRBFGenerator(dimension, partials, f_func)
end

function get_field{T}(gen::HRBFGenerator, points::Array{T}, normals::Array{T})
	dimension = gen.dimension
	@assert dimension == size(points, 2)
	@assert dimension == size(normals, 2)
	num_points = size(points, 1)

	A_type = promote_type(T, Float64)
	A_11 = A_type[gen.partials[(:f, :lambda)](i,j)(points) for i in 1:num_points, j in 1:num_points]
	A_12 = hcat([A_type[gen.partials[(:f, :beta)](i,j,k)(points) for i in 1:num_points, k in 1:dimension] for j in 1:num_points]...)
	A_21 = vcat([A_type[gen.partials[(:g, :lambda)](i, j, k)(points) for k in 1:dimension, j in 1:num_points] for i in 1:num_points]...)
	A_22 = vcat([hcat([A_type[gen.partials[(:g, :beta)](i, j, k, l)(points) for k in 1:dimension, l in 1:dimension] for j in 1:num_points]...) for i in 1:num_points]...)

	A = [A_11 A_12; A_21 A_22]
	A[map(isnan, A)] = 0.0
	b = vcat(zeros(num_points), normals'[:])
	y = A \ b
	lambda_values = y[1:num_points]
	beta_values = reshape(y[num_points+1:end], dimension, num_points)'
	return (x_values...) -> sum([gen.field_func(lambda_values[i], beta_values[i,:]..., x_values..., points[i,:]...) for i = 1:num_points])
end

end # module


text:
using CSV
using DataFrames
using Plots
using EcologicalNetworks
using EcologicalNetworksPlots
using Distributions

function load_data()
    df = CSV.read("WoodEtal_Append1_v2.csv")

    # Separate based on WebScale
    scale = [df[df.WebScale .== u,:] for u in unique(df.WebScale)]

    # select for edge interaction columns only
    edgelist_T = scale[1][:, [:PredTSN, :PreyTSN]]

    # view all entries that contain "san" in PredTSN
    pred_san = edgelist_T[startswith.(edgelist_T.PredTSN, "san"),:]
    # remove "san"
    pred_san.PredTSN = replace.(pred_san.PredTSN, "san" .=> "")
    pred_san
    #view unique predID without san
    unique(pred_san,1)
    edgelist_T.PredTSN == "102" #manually check if in original df

    # remove the "san" from the dataframe itself
     edgelist_T.PredTSN = replace.(edgelist_T.PredTSN, "san" .=> "")

    # Do the same for PreyTSN
    backup2 = copy(edgelist_T)
    old_ids2 = unique(backup2.PreyTSN)
    # keep the indexes before the change
    prey_indexes = startswith.(edgelist_T.PreyTSN, "san")
    # remove the "san" from the dataframe itself
    edgelist_T.PreyTSN = replace.(edgelist_T.PreyTSN, "san" .=> "")
    edgelist_T
    #check unique preyIDs without san
    check_uniques2 = [in(new_id2, old_ids2) for new_id2 in edgelist_T.PreyTSN[prey_indexes]]
    filter(isone, check_uniques2)

    # convert "string" to "integer"
    edgelist_T.PredTSN = parse.(Int64, edgelist_T.PredTSN)
    edgelist_T.PreyTSN = parse.(Int64, edgelist_T.PreyTSN)

    return edgelist_T
end

max_stars_repo_name: kiristern/reseaux

max_issues_repo_name: kiristern/reseaux

max_forks_repo_name: kiristern/reseaux


text:
# Formulas for representing and working with linear-model-type expressions
# Original by Harlan D. Harris.  Later modifications by John Myles White
# and Douglas M. Bates.

## Formulas are written as expressions and parsed by the Julia parser.
## For example :(y ~ a + b + log(c))
## In Julia the & operator is used for an interaction.  What would be written
## in R as y ~ a + b + a:b is written :(y ~ a + b + a&b) in Julia.
## The equivalent R expression, y ~ a*b, is the same in Julia

## The lhs of a one-sided formula is 'nothing'
## The rhs of a formula can be 1

type Formula
    lhs::@compat(Union{Symbol, Expr, Void})
    rhs::@compat(Union{Symbol, Expr, Integer})
end

macro ~(lhs, rhs)
    ex = Expr(:call,
              :Formula,
              Base.Meta.quot(lhs),
              Base.Meta.quot(rhs))
    return ex
end

type Terms
    terms::Vector
    eterms::Vector        # evaluation terms
    factors::Matrix{Int8} # maps terms to evaluation terms
    order::Vector{Int}    # orders of rhs terms
    response::Bool        # indicator of a response, which is eterms[1] if present
    intercept::Bool       # is there an intercept column in the model matrix?
end

type ModelFrame
    df::AbstractDataFrame
    terms::Terms
    msng::BitArray
end

type ModelMatrix{T <: @compat(Union{Float32, Float64})}
    m::Matrix{T}
    assign::Vector{Int}
end

Base.size(mm::ModelMatrix) = size(mm.m)
Base.size(mm::ModelMatrix, dim...) = size(mm.m, dim...)

function Base.show(io::IO, f::Formula)
    print(io,
          string("Formula: ",
                 f.lhs == nothing ? "" : f.lhs, " ~ ", f.rhs))
end

## Return, as a vector of symbols, the names of all the variables in
## an expression or a formula
function allvars(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    cc=Symbol[]
    for i in ex.args[2:end] cc=append!(cc,allvars(i)) end
    cc
end
allvars(f::Formula) = unique(vcat(allvars(f.rhs), allvars(f.lhs)))
allvars(sym::Symbol) = [sym]
allvars(v::Any) = Array(Symbol, 0)

# special operators in formulas
const specials = Set([:+, :-, :*, :/, :&, :|, :^])

function dospecials(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    a1 = ex.args[1]
    if !(a1 in specials) return ex end
    excp = copy(ex)
    excp.args = vcat(a1,map(dospecials, ex.args[2:end]))
    if a1 != :* return excp end
    aa = excp.args
    a2 = aa[2]
    a3 = aa[3]
    if length(aa) > 3
        excp.args = vcat(a1, aa[3:end])
        a3 = dospecials(excp)
    end
    ## this order of expansion gives the R-style ordering of interaction
    ## terms (after sorting in increasing interaction order) for higher-
    ## order interaction terms (e.g. x1 * x2 * x3 should expand to x1 +
    ## x2 + x3 + x1&x2 + x1&x3 + x2&x3 + x1&x2&x3)
    :($a2 + $a2 & $a3 + $a3)
end
dospecials(a::Any) = a

## Distribution of & over +
const distributive = @compat Dict(:& => :+)

distribute(ex::Expr) = distribute!(copy(ex))
distribute(a::Any) = a
## apply distributive property in-place
function distribute!(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    [distribute!(a) for a in ex.args[2:end]]
    ## check that top-level can be distributed
    a1 = ex.args[1]
    if a1 in keys(distributive)

        ## which op is being DISTRIBUTED (e.g. &, *)?
        distributed_op = a1
        ## which op is doing the distributing (e.g. +)?
        distributing_op = distributive[a1]

        ## detect distributing sub-expression (first arg is, e.g. :+)
        is_distributing_subex(e) =
            typeof(e)==Expr && e.head == :call && e.args[1] == distributing_op

        ## find first distributing subex
        first_distributing_subex = findfirst(is_distributing_subex, ex.args)
        if first_distributing_subex != 0
            ## remove distributing subexpression from args
            subex = splice!(ex.args, first_distributing_subex)

            newargs = Any[distributing_op]
            ## generate one new sub-expression, which calls the distributed operation
            ## (e.g. &) on each of the distributing sub-expression's arguments, plus
            ## the non-distributed arguments of the original expression.
            for a in subex.args[2:end]
                new_subex = copy(ex)
                push!(new_subex.args, a)
                ## need to recurse here, in case there are any other
                ## distributing operations in the sub expression
                distribute!(new_subex)
                push!(newargs, new_subex)
            end
            ex.args = newargs
        end
    end
    ex
end
distribute!(a::Any) = a

const associative = Set([:+,:*,:&])       # associative special operators

## If the expression is a call to the function s return its arguments
## Otherwise return the expression
function ex_or_args(ex::Expr,s::Symbol)
    if ex.head != :call error("Non-call expression encountered") end
    if ex.args[1] == s
        ## recurse in case there are more :calls of s below
        return vcat(map(x -> ex_or_args(x, s), ex.args[2:end])...)
    else
        ## not a :call to s, return condensed version of ex
        return condense(ex)
    end
end
ex_or_args(a,s::Symbol) = a

## Condense calls like :(+(a,+(b,c))) to :(+(a,b,c))
function condense(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    a1 = ex.args[1]
    if !(a1 in associative) return ex end
    excp = copy(ex)
    excp.args = vcat(a1, map(x->ex_or_args(x,a1), ex.args[2:end])...)
    excp
end
condense(a::Any) = a

## always return an ARRAY of terms
getterms(ex::Expr) = (ex.head == :call && ex.args[1] == :+) ? ex.args[2:end] : Expr[ex]
getterms(a::Any) = Any[a]

ord(ex::Expr) = (ex.head == :call && ex.args[1] == :&) ? length(ex.args)-1 : 1
ord(a::Any) = 1

const nonevaluation = Set([:&,:|])        # operators constructed from other evaluations
## evaluation terms - the (filtered) arguments for :& and :|, otherwise the term itself
function evt(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    if !(ex.args[1] in nonevaluation) return ex end
    filter(x->!isa(x,Number), vcat(map(getterms, ex.args[2:end])...))
end
evt(a) = Any[a]

function Terms(f::Formula)
    rhs = condense(distribute(dospecials(f.rhs)))
    tt = unique(getterms(rhs))
    tt = tt[!(tt .== 1)]             # drop any explicit 1's
    noint = (tt .== 0) | (tt .== -1) # should also handle :(-(expr,1))
    tt = tt[!noint]
    oo = Int[ord(t) for t in tt]     # orders of interaction terms
    if !issorted(oo)                 # sort terms by increasing order
        pp = sortperm(oo)
        tt = tt[pp]
        oo = oo[pp]
    end
    etrms = map(evt, tt)
    haslhs = f.lhs != nothing
    if haslhs
        unshift!(etrms, Any[f.lhs])
        unshift!(oo, 1)
    end
    ev = unique(vcat(etrms...))
    sets = [Set(x) for x in etrms]
    facs = Int8[t in s for t in ev, s in sets]
    Terms(tt, ev, facs, oo, haslhs, !any(noint))
end

function remove_response(t::Terms)
    # shallow copy original terms
    t = Terms(t.terms, t.eterms, t.factors, t.order, t.response, t.intercept)
    if t.response
        t.order = t.order[2:end]
        t.eterms = t.eterms[2:end]
        t.factors = t.factors[2:end, 2:end]
        t.response = false
    end
    return t
end

## Default NA handler.  Others can be added as keyword arguments
function na_omit(df::DataFrame)
    cc = complete_cases(df)
    df[cc,:], cc
end

## Trim the pool field of da to only those levels that occur in the refs
function dropUnusedLevels!(da::PooledDataArray)
    rr = da.refs
    uu = unique(rr)
    length(uu) == length(da.pool) && return da
    T = eltype(rr)
    su = sort!(uu)
    dict = Dict(zip(su, one(T):convert(T, length(uu))))
    da.refs = map(x -> dict[x], rr)
    da.pool = da.pool[uu]
    da
end
dropUnusedLevels!(x) = x

function ModelFrame(trms::Terms, d::AbstractDataFrame)
    df, msng = na_omit(DataFrame(map(x -> d[x], trms.eterms)))
    names!(df, convert(Vector{Symbol}, map(string, trms.eterms)))
    for c in eachcol(df) dropUnusedLevels!(c[2]) end
    ModelFrame(df, trms, msng)
end

ModelFrame(f::Formula, d::AbstractDataFrame) = ModelFrame(Terms(f), d)
ModelFrame(ex::Expr, d::AbstractDataFrame) = ModelFrame(Formula(ex), d)

function StatsBase.model_response(mf::ModelFrame)
    mf.terms.response || error("Model formula one-sided")
    convert(Array, mf.df[round(Bool, mf.terms.factors[:, 1])][:, 1])
end

function contr_treatment(n::Integer, contrasts::Bool, sparse::Bool, base::Integer)
    if n < 2 error("not enought degrees of freedom to define contrasts") end
    contr = sparse ? speye(n) : eye(n) .== 1.
    if !contrasts return contr end
    if !(1 <= base <= n) error("base = $base is not allowed for n = $n") end
    contr[:,vcat(1:(base-1),(base+1):end)]
end
contr_treatment(n::Integer,contrasts::Bool,sparse::Bool) = contr_treatment(n,contrasts,sparse,1)
contr_treatment(n::Integer,contrasts::Bool) = contr_treatment(n,contrasts,false,1)
contr_treatment(n::Integer) = contr_treatment(n,true,false,1)
cols(v::PooledDataVector) = contr_treatment(length(v.pool))[v.refs,:]
cols(v::DataVector) = convert(Vector{Float64}, v.data)
cols(v::Vector) = convert(Vector{Float64}, v)

function isfe(ex::Expr)                 # true for fixed-effects terms
    if ex.head != :call error("Non-call expression encountered") end
    ex.args[1] != :|
end
isfe(a) = true

function expandcols(trm::Vector)
    if length(trm) == 1
        return convert(Array{Float64}, trm[1])
    else
        a = convert(Array{Float64}, trm[1])
        b = expandcols(trm[2:end])
        nca = size(a, 2)
        ncb = size(b, 2)
        return hcat([a[:, i] .* b[:, j] for i in 1:nca, j in 1:ncb]...)
    end
end

function nc(trm::Vector)
    isempty(trm) && return 0
    n = 1
    for x in trm
        n *= size(x, 2)
    end
    n
end

function ModelMatrix(mf::ModelFrame)
    trms = mf.terms
    aa = Any[Any[ones(size(mf.df,1), @compat(Int(trms.intercept)))]]
    asgn = zeros(Int, @compat(Int(trms.intercept)))
    fetrms = Bool[isfe(t) for t in trms.terms]
    if trms.response unshift!(fetrms, false) end
    ff = trms.factors[:, fetrms]
    ## need to be cautious here to avoid evaluating cols for a factor with many levels
    ## if the factor doesn't occur in the fetrms
    rows = Bool[x != 0 for x in sum(ff, 2)]
    ff = ff[rows, :]
    cc = [cols(col) for col in columns(mf.df[:, rows])]
    for j in 1:size(ff,2)
        trm = cc[round(Bool, ff[:, j])]
        push!(aa, trm)
        asgn = vcat(asgn, fill(j, nc(trm)))
    end
    ModelMatrix{Float64}(hcat([expandcols(t) for t in aa]...), asgn)
end

termnames(term::Symbol, col) = [string(term)]
function termnames(term::Symbol, col::PooledDataArray)
    levs = levels(col)
    [string(term, " - ", levs[i]) for i in 2:length(levs)]
end

function coefnames(mf::ModelFrame)
    if mf.terms.intercept
        vnames = UTF8String["(Intercept)"]
    else
        vnames = UTF8String[]
    end
    # Need to only include active levels
    for term in mf.terms.terms
        if isa(term, Expr)
            if term.head == :call && term.args[1] == :|
                continue                # skip random-effects terms
            elseif term.head == :call && term.args[1] == :&
                ## for an interaction term, combine term names pairwise,
                ## starting with rightmost terms
                append!(vnames,
                        foldr((a,b) ->
                              vec([string(lev1, " & ", lev2) for
                                   lev1 in a,
                                   lev2 in b]),
                              map(x -> termnames(x, mf.df[x]), term.args[2:end])))
            else
                error("unrecognized term $term")
            end
        else
            append!(vnames, termnames(term, mf.df[term]))
        end
    end
    return vnames
end

function Formula(t::Terms)
    lhs = t.response ? t.eterms[1] : nothing
    rhs = Expr(:call,:+)
    if t.intercept
        push!(rhs.args,1)
    end
    append!(rhs.args,t.terms)
    Formula(lhs,rhs)
end


text:
using Gridap
using Gridap.Io
using GridapGmsh

model = GmshDiscreteModel("elasticFlag_coarse.msh")

writevtk(model,"elasticFlag_coarse")

fn = "elasticFlag_coarse.json"
to_json_file(model,fn)


text:
module PkgTest

include("df_test.jl")
include("normal_normal_lognormal_test.jl")
include("mvt_gamma_invgamma_test.jl")

# include("mvcalc_test.jl")
# include("perf_test.jl")
# include("startingpoint_test.jl")

end # PkgTest module


text:
SUPPORTED_REALS = (
    Float32, Float64,
    Int8, Int16, Int32, Int64,
    UInt8, UInt16, UInt32, UInt64,
    Bool,
    )
SUPPORTED_COMPLEX = (ComplexF32, ComplexF64)

function int2bytar(i::Int; i_type = UInt32)
    if i_type == Bool
        return UInt8(Bool)
    else
        return reinterpret(UInt8, [i_type(i)])
    end
end

function bytar2int(b::ByteArr; i_type = UInt32)
    if i_type != Bool
        return reinterpret(i_type, b)[1]
    else
        return Bool(b[1])
    end
end

function numtypestring(ar)
    t = eltype(ar)
    realtypes =
        (Float32, Float64, Int8, Int16, Int32, Int64, UInt8, UInt16, UInt32, UInt64, Bool)
    if t in realtypes
        return string(t)
    elseif t == ComplexF32  # special cases, as string(ComplexF32) is "Complex{Float32}"
        return "ComplexF32"
    elseif t == ComplexF64
        return "ComplexF64"
    else
        return throw(DomainError("$(string(t)) is not a supported numeric type for exchange with LabVIEW: consider converting array first."))
    end
end

function fromrowmajor(v, arrdims)
    revdims = Tuple(reverse(arrdims))
    neworder = length(revdims):-1:1
    arr = permutedims(reshape(v, revdims), neworder)
    return arr
end

function fromcolmajor(v, arrdims)
    revdims = Tuple(reverse(arrdims))
    neworder = length(revdims):-1:1
    arr = permutedims(reshape(v, revdims), neworder)
    return arr
end

function cmplxswap(c::Complex)
    c = imag(c) + real(c)im
    return c
end

function emptyarr(arrdims, numtype)
    Array{numtype}(undef, zeros(Int, length(arrdims))...)
end

function bin2num(; bin_data, nofbytes, start, arrdims, numtype)
    nt = Symbol(numtype)
    numtype = eval(nt)

    if nofbytes == 0 || any(arrdims.==0)
        return emptyarr(arrdims, numtype)
    end

    bin_data = bin_data[start:start+nofbytes-1]

    if numtype != Bool
        nums = collect(reinterpret(numtype, bin_data))
    else
        nums = Bool.(bin_data)
    end

    if length(arrdims) > 1
        nums = fromrowmajor(nums, arrdims)
    end

    if eltype(nums) in (ComplexF32, ComplexF64)
        nums .= cmplxswap.(nums)
    end
    return nums
end

function isimage(a)
    return (typeof(a) <: AbstractArray{C,2} where C <: Color)
end

function bin2img(; bin_data, nofbytes, start, arrdims, numtype)

    if nofbytes == 0 || any(arrdims.==0)
        T = typeof(RGB(0,0,0)) # type by example
        return emptyarr(arrdims, T)
    end

    bin_data = bin_data[start:start+nofbytes-1]

    if numtype == "img24bit"
        arrdims = (3, arrdims[2], arrdims[1])
        nums = bin_data
        nums = reshape(bin_data, arrdims)
        nums = permutedims(nums, [1, 3, 2])
        nums = colorview(RGB, normedview(nums))
    else
        throw(DomainError("image data type ($numtype) not supported"))
    end

    return nums
end

function fn_by_category(ctg)
    fs = Dict("images" => bin2img, "numarrays" => bin2num)
    return fs[ctg]
end

"""
   bin2nums(; bin_data, bindata_descr)

Get data as received from LabVIEW (binary + array of decriptions) and return a
dictionary of de-serialized arrays (or other supported data types, e.g. images).
The dictionary items meant to be supplied as kwargs to user functions.
"""
function bin2nums(; bin_data, bindata_descr)
    arrs = [
        fn_by_category(bdd.category)(
            bin_data = bin_data,
            nofbytes = bdd.nofbytes,
            start = bdd.start,
            arrdims = bdd.arrdims,
            numtype = bdd.numtype,
        ) for bdd in bindata_descr
    ]
    kwarg_names = [Symbol(bdd.kwarg_name) for bdd in bindata_descr]
    darrs = Dict(Pair.(kwarg_names, arrs))
    return darrs
end

# # # # # # # #

function tocolmajor(arr)
    arrdims = size(arr)
    lng = length(arrdims)
    if lng > 1
        neworder = lng:-1:1
        arr = permutedims(arr, neworder)
    end
    return vec(arr)
end

function nums2ByteArr(nums)
    nums = tocolmajor(nums)
    if eltype(nums) == Bool
        return UInt8.(nums)
    else
        return collect(reinterpret(UInt8, nums))
    end
end

function data2bin(img::I, kwarg_name) where I <: AbstractArray{C,2} where C <: Colorant

    # in most cases, image is already RGB{N0f8} anyway, but who knows:
    img = convert.(RGB{N0f8}, img)

    cnv = channelview(img)
    rwv = collect(rawview(cnv))
    prw = collect(permutedims(rwv, [1,3,2]))
    bd = vec(prw)

    bdd = BinDescr()
    bdd.kwarg_name = kwarg_name
    bdd.numtype = "img24bit"
    bdd.category = "images"
    bdd.arrdims = collect(size(img))
    bdd.nofbytes = length(bd)
    return (bd, bdd)
end

function data2bin(arrdata, kwarg_name)
    bdd = BinDescr()
    bdd.kwarg_name = kwarg_name
    bdd.numtype = numtypestring(arrdata)
    if eltype(arrdata) <: Complex
        arrdata .= cmplxswap.(arrdata)
    end
    bdd.arrdims = collect(size(arrdata))
    bd = nums2ByteArr(arrdata)
    bdd.nofbytes = length(bd)
    return (bd, bdd)
end

"""
   data2bin!!(bin_data::ByteArr, bindata_descr::Vector{BinDescr}; arrdata, kwarg_name)

Serialize for exchange with LabVIEW: Compute binary representation and
description of data (currently supported types are numeric arrays and images).
Append the binary representation to `bin_data`, and description to
`bindata_descr` array (if the first two arguments are provided, mutate them,
otherwise create).
"""
function data2bin!!(
    bin_data::ByteArr = ByteArr(),
    bindata_descr::Vector{BinDescr} = BinDescr[];
    arrdata,
    kwarg_name,
)
    @assert isempty(bin_data) == isempty(bindata_descr)
    bd, bdd = data2bin(arrdata, kwarg_name)
    bdd.start = length(bin_data) + 1
    bin_data = vcat(bin_data, bd)
    push!(bindata_descr, bdd)
    return (; bin_data, bindata_descr)
end

"""
   nums2bin(arrs)

Accept an iterable of array-like data (numeric arrays, images, anything to come...). Compute
binary representation and description for them and return binary data and array of
descriptions, ready to be sent to LabVIEW.
"""
function nums2bin(arrs)
    bin_data = ByteArr()
    bindata_descr = BinDescr[]

    for arrname in keys(arrs)
        bin_data, bindata_descr =
            data2bin!!(bin_data, bindata_descr; arrdata = arrs[arrname], kwarg_name = arrname)
    end
    return (; bin_data, bindata_descr)
end


text:
"""
    weibull_rng(λ, β, shape=1; seed=nothing)

Generate a `shape` element array of random variables from a Weibull(`λ`, `β`) distribution. Optionally you can set a specific seed.

# Notes

The pdf of an Weibull(λ, β) distribution is given as:

``f(x, λ, β) = λ e^{-λx^β} \\quad x ≥ 0``

# Examples

```julia-repl
julia> weibull_rng(2, 2, seed=42)
1-element Vector{Float64}:
 0.31445725834527055

julia> weibull_rng(2, 2, 2)
2-element Vector{Float64}:
 0.39561285703154575
 0.6021921673483441

julia> weibull_rng(2, 2, (2,2))
2×2 Matrix{Float64}:
 0.428896  0.109897
 0.812854  0.427906

```

# References

C. Alexopoulos, D. Goldsman. Random variate generation. 2020.

Law, A. Simulation modeling and analysis, 5th Ed. McGraw Hill Education, Tuscon, 2013.
"""
function weibull_rng(λ::Real, β::Real, shape::Union{Int, Tuple{Vararg{Int}}}=1; seed::Union{Int, Nothing}=nothing)
    U = get_std_uniform(shape, seed=seed)
    X = (1/λ) .* (-log.(1 .- U)) .^ (1/β)
    return X
end


text:
""
function run_opf_variable_impedance_loss(file, optimizer; kwargs...)
    return run_opf_variable_impedance_loss(file, _PM.ACPPowerModel, optimizer; kwargs...)
end

""
function run_opf_variable_impedance_loss(file, model_type::Type, optimizer; kwargs...)
    return _PM.run_model(file, model_type, optimizer, build_opf_variable_impedance_loss; kwargs...)
end

"This is our new model in accordance with Algorithm 1 of the paper"
function build_opf_variable_impedance_loss(pm::_PM.AbstractPowerModel)
    # Comments reference slide 42 on installfest pdf
    _PM.variable_bus_voltage(pm) # Variable in V_i in model 1 -> Constraint 1.3
    _PM.variable_gen_power(pm) # Variable S_i in model 1 -> Constraint 1.5
    _PM.variable_branch_power(pm) # Implicit
    # _PM.variable_dcline_flow(pm)
    variable_admittance(pm) # s4, s5 set upper and lower bounds of g, b

    minimum_impedance_distance(pm) # Our new objective

    constraint_loss_faithfulness(pm) # s2
    _PM.constraint_model_voltage(pm) # Implicit

    for i in _PM.ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i) # Implicit
    end

    for i in _PM.ids(pm, :bus)
        _PM.constraint_power_balance(pm, i) # Does this reference constraint_kcl_shunt ie Constraint 1.7?
    end

    for i in _PM.ids(pm, :branch)
        constraint_ohms_from_variable_impedance(pm, i) # Constraint s3, s4, 1.8
        constraint_ohms_to_variable_impedance(pm, i) # Constraint s3, s4, 1.8

        _PM.constraint_voltage_angle_difference(pm, i) # Constraint 1.4

        _PM.constraint_thermal_limit_from(pm, i) # Constraint 1.6
        _PM.constraint_thermal_limit_to(pm, i) # Constraint 1.6
    end

    # for i in _PM.ids(pm, :dcline)
    #     _PM.constraint_dcline(pm, i)
    # end
end

""
function run_opf_variable_impedance_cost(file, optimizer; kwargs...)
    return run_opf_variable_impedance_cost(file, _PM.ACPPowerModel, optimizer; kwargs...)
end

""
function run_opf_variable_impedance_cost(file, model_type::Type, optimizer; kwargs...)
    return _PM.run_model(file, model_type, optimizer, build_opf_variable_impedance_cost; kwargs...)
end

"This is our new model in accordance with Algorithm 1 of the paper"
function build_opf_variable_impedance_cost(pm::_PM.AbstractPowerModel)
    # Comments reference slide 42 on installfest pdf
    _PM.variable_bus_voltage(pm) # Variable in V_i in model 1 -> Constraint 1.3
    _PM.variable_gen_power(pm) # Variable S_i in model 1 -> Constraint 1.5
    _PM.variable_branch_power(pm) # Implicit
    # _PM.variable_dcline_flow(pm)
    variable_admittance(pm) # s4, s5 set upper and lower bounds of g, b
    variable_fuel_cost(pm)

    minimum_impedance_distance(pm) # Our new objective

    constraint_fuel_cost_quadratic(pm)
    constraint_cost_faithfulness(pm) # s2

    _PM.constraint_model_voltage(pm) # Implicit

    for i in _PM.ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i) # Implicit
    end

    for i in _PM.ids(pm, :bus)
        _PM.constraint_power_balance(pm, i) # Does this reference constraint_kcl_shunt ie Constraint 1.7?
    end

    for i in _PM.ids(pm, :branch)
        constraint_ohms_from_variable_impedance(pm, i) # Constraint s3, s4, 1.8
        constraint_ohms_to_variable_impedance(pm, i) # Constraint s3, s4, 1.8

        _PM.constraint_voltage_angle_difference(pm, i) # Constraint 1.4

        _PM.constraint_thermal_limit_from(pm, i) # Constraint 1.6
        _PM.constraint_thermal_limit_to(pm, i) # Constraint 1.6
    end

    # for i in _PM.ids(pm, :dcline)
    #     _PM.constraint_dcline(pm, i)
    # end
end


text:
using Markdown
using InteractiveUtils

# ╔═╡ f9fb3324-bf08-11eb-3893-85f5df3d567d
using Revise

# ╔═╡ 216ccd98-4953-4d9d-b770-69caee1c3fcc
using Reproduce, ReproducePlotUtils, StatsPlots, RollingFunctions, Statistics, FileIO, PlutoUI, Pluto

# ╔═╡ 50c823af-0573-4be3-98b8-b40cc1c8e501
color_scheme = [
    colorant"#44AA99",
    colorant"#332288",
    colorant"#DDCC77",
    colorant"#999933",
    colorant"#CC6677",
    colorant"#AA4499",
    colorant"#DDDDDD",
	colorant"#117733",
	colorant"#882255",
	colorant"#88CCEE",
]

# ╔═╡ 44b42144-1ca4-4814-993f-6984687e3c7b
const RPU = ReproducePlotUtils

# ╔═╡ ded0bd85-55b5-4d22-8ad2-b510381aadae
ic_sr, dd_sr = RPU.load_data("../data/ER/OneDTMaze_RR_ER_cor_err/")

# ╔═╡ 303bfc74-2abc-492f-89ef-c5f15322cfec
data_sr_avg = RPU.get_line_data_for(
		ic_sr,
        ["batch_size"],
        ["alpha_init", "eta"];
        comp=:min,
		get_comp_data=(x)->sum(x["results"][:oned_tmaze_dmu_error][:, end-30]),
        get_data=(x)->mean(x["results"][:oned_tmaze_dmu_error], dims=1)[1, :])

# ╔═╡ e9527191-dcac-43b6-aafc-d8a7f9750917
data_sr = [RPU.get_line_data_for(
		ic_sr,
        ["batch_size"],
        ["alpha_init", "eta"];
        comp=:min,
		get_comp_data=(x)->sum(x["results"][:oned_tmaze_dmu_error][:, end-30]),
        get_data=(x)->x["results"][:oned_tmaze_dmu_error][i,:]) for i in 1:4]

# ╔═╡ d508ab1b-8016-459e-ba70-7553ca26ae04
let
	tls = ["Distractor", "Constant-1", "Drifter", "Constant-2"]
	plts = [plot(data_sr[i], label_idx="batch_size",
		 legendtitle=nothing, palette=color_scheme, title=tls[i], size=(800,800)) for i in 1:4]
	plot(plts...)
end

# ╔═╡ db7ffcb4-e6f9-4591-a35a-9fb4b7d5925f
FileIO.load(joinpath(ic_sr[1].folder_str, "results.jld2"))

# ╔═╡ c427e9d6-dbe7-49ae-bfc9-164612233a5d
ic_q, dd_q = RPU.load_data("../local_data/OneDTMaze_RR_ER_2/")

# ╔═╡ 151330f1-b8b5-4a62-be6d-37fe0870e263
data_q_avg = RPU.get_line_data_for(
		ic_q,
        ["batch_size"],
        ["alpha_init", "eta"];
        comp=:min,
		get_comp_data=(x)->sum(x["results"][:oned_tmaze_dmu_error][:, end-30]),
        get_data=(x)->mean(x["results"][:oned_tmaze_dmu_error], dims=1)[1, :])

# ╔═╡ d06a43c2-b4e7-48c7-b844-4a12ce7b8d98
let
	bs = 1
	plot(data_sr_avg, Dict("batch_size"=>bs), label="SR $(bs)", label_idx="batch_size",
		 legendtitle=nothing, palette=color_scheme)
	plot!(data_q_avg, Dict("batch_size"=>bs), label="Q $(bs)", label_idx="batch_size",
		 legendtitle=nothing, palette=color_scheme)
end

# ╔═╡ 0744dfcf-d695-4257-a2d1-a2433cc96428
data_q = [RPU.get_line_data_for(
		ic_q,
        ["batch_size"],
        ["alpha_init", "eta"];
        comp=:min,
		get_comp_data=(x)->sum(x["results"][:oned_tmaze_dmu_error][:, end-30]),
        get_data=(x)->x["results"][:oned_tmaze_dmu_error][i,:]) for i in 1:4]

# ╔═╡ 4e168ab6-4cc0-48c5-a5f4-ade5e4b87746
let
	tls = ["Distractor", "Constant-1", "Drifter", "Constant-2"]
	plts = [plot(data_q[i], label_idx="batch_size",
		 legendtitle=nothing, palette=color_scheme, title=tls[i], size=(800,800)) for i in 1:4]
	plot(plts...)
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
module LSystem

export @lsys
export LModel, add_rule!
export LState, next, result

using MacroTools

# L-System model definition

"""
A L-system model is represented by an axiom called `axiom`
and a set of rewriting `rules`.
"""
struct LModel
    axiom
    rules
end

"Create a L-system model."
LModel(axiom) = LModel([axiom], Dict())

"Add rule to a model."
function add_rule!(model::LModel, left::T, right::T) where {T <: AbstractString}
    model.rules[left] = split(right, "")
    return nothing
end

"Display model nicely."
function Base.show(io::IO, model::LModel)
    println(io, "LModel:")
    println(io, "  Axiom: ", join(model.axiom))
    for k in sort(collect(keys(model.rules)))
        println(io, "  Rule:  ", k, " → ", join(model.rules[k]))
    end
end

# Tracking state of the system

"""
A L-system state contains a reference to the `model`, the
current iteration, and the result.
"""
struct LState
    model
    current_iteration
    result
end

"Create a L-system state from a `model`."
LState(model::LModel) = LState(model, 1, model.axiom)

"Advance to the next state and returns a new LState object."
function next(state::LState)
    new_result = []
    for el in state.result
        next_el = get(state.model.rules, el, el)
        push!.(Ref(new_result), next_el)
    end
    return LState(state.model, state.current_iteration + 1, new_result)
end

"Repeated next call"
next(state::LState, n) = n > 0 ? next(next(state), n-1) : state

"Compact the result suitable for display"
result(state::LState) = join(state.result)

Base.show(io::IO, s::LState) =
    print(io, "LState(", s.current_iteration, "): ", result(s))

# DSL implementation

"""
The @lsys macro is used to construct a L-System model object [LModel](@ref).
The domain specific language requires a single axiom and a set of rewriting rules.
For example:

```
model = @lsys begin
    axiom : A
    rule  : A → AB
    rule  : B → A
end
```
"""
macro lsys(ex)
    ex = MacroTools.postwalk(walk, ex)
    push!(ex.args, :( model ))
    return ex
end

# Walk the AST tree and match expressions.
function walk(ex)

    match_axiom = @capture(ex, axiom : sym_)
    if match_axiom
        sym_str = String(sym)
        return :( model = LModel($sym_str) )
    end

    match_rule = @capture(ex, rule : original_ → replacement_)
    if match_rule
        original_str = String(original)
        replacement_str = String(replacement)
        return :(
            add_rule!(model, $original_str, $replacement_str)
        )
    end

    return ex
end

end # module


text:
# thsi script calculates the barycentric coordinates for 2D and 3D
# SBP operator nodes, and prints the case statement needed by
# apf::FieldShape

using SummationByParts
include("nodecalc3.jl")

function nodecalc(sbp::TriSBP, isDG::Bool)
  vtx = [0.0 0; 1 0; 0 1]

  if isDG
    coords = SummationByParts.SymCubatures.calcnodes(sbp.cub, vtx)
  else
    coords = calcnodes(sbp, vtx)
  end

  xi = convertToBary(coords)

  return xi, coords
end

function convertToBary(coords_xy::Array)
# convert x-y coordinates to barycentric coordinates
  dim = size(coords_xy, 1)

  if dim == 2
    vtx = [0.0 0; 1 0; 0 1]
    r1 = vtx[1, :]
    r2 = vtx[2, :]
    r3 = vtx[3, :]

    T = zeros(2,2)
    T[:, 1] = r2 - r1
    T[:, 2] = r3 - r1
  elseif dim == 3

    vtx = [-1.0 -1.0 -1.0
          1.0 -1.0 -1.0
         -1.0 1.0 -1.0
         -1.0 -1.0 1.0]
    #=
    vtx = [0 0  0;
           1. 0 0;
           0 1 0;
           0 0 1]
    =#
    r1 = vtx[1, :]
    r2 = vtx[2, :]
    r3 = vtx[3, :]
    r4 = vtx[4, :]
    T = zeros(3,3)
    T[:, 1] = r2 - r1
    T[:, 2] = r3 - r1
    T[:, 3] = r4 - r1
    println("T = \n", T)
    println("r1 = \n", r1)
    println("size(T) = ", size(T))
    println("size(r1 = ", size(r1))
  end

  xi = zeros(coords_xy)
  for i=1:size(coords_xy,2)
    xi[:, i] = T\(coords_xy[:, i] - r1)
  end

  return xi
end

function nodecalc(sbp::TetSBP, isDG::Bool)

  vtx = [-1.0 -1.0 -1.0
        1.0 -1.0 -1.0
       -1.0 1.0 -1.0
       -1.0 -1.0 1.0]

  if isDG
    coords = SummationByParts.SymCubatures.calcnodes(sbp.cub, vtx)
  else
    coords = calcnodes(sbp, vtx)
  end

  xi = convertToBary(coords)

  return xi, coords
end

function minNodeDist(sbp, isDG::Bool)
# get the minimum distance between nodes on a reference element of degree p

    xi, coords = nodecalc(sbp, isDG)
    min_dist = typemax(Float64)
#    println("coords = ", coords)
    for i=1:size(coords, 2)
      for j=(i+1):size(coords, 2)
	# calculate distance between node i and node j
	dist_j = norm(coords[:, i] - coords[:, j])

	if dist_j < min_dist
	  min_dist = dist_j
	end

      end  # end loop j
    end  # end loop i

#    println("for p=$p elements, min node distance = ", min_dist)

  return min_dist
end

#sbp = getTriSBPWithDiagE(degree=4, vertices=false)
#sbp = getTetSBPWithDiagE(degree=2)a
#sbp = getTriSBPOmega2(degree=4)
#sbp = TriSBP{Float64}(degree=1, reorder=false, internal=false)
#sbp = TetSBP{Float64}(degree=3, reorder=false, internal=true)
#sbp = getTetSBPDiagE(degree=4)
sbp = getTriSBPDiagE(degree=4, mincond=false)
#sbp = getTriSBPOmega(degree=3)
xi, coords = nodecalc(sbp, true)
printCaseStatement(xi)
printBaryCoords(xi, coords)
#writedlm("coordsout.dat", coords, ' ')
#minNodeDist(2)

#=
coords = [0.5 0 0;
          0 1 0;
          0 0 1]
xi = convertToBary(coords)
println("xi = \n", xi)
=#


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

using Markdown
using InteractiveUtils

# ╔═╡ 1c31a6c8-fe76-11ea-1122-0580af80a98b
using Pkg, DrWatson

# ╔═╡ 1c31e480-fe76-11ea-356d-6b2d21e37ace
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanQuap
	using StatisticalRethinking
end

# ╔═╡ d1b31aa0-fe75-11ea-0d96-bdbf19c2c1c7
md"## Clip-06-06-09s.jl"

# ╔═╡ 1c325ef6-fe76-11ea-1967-45317ad16daa
begin
	df = CSV.read(sr_datadir("milk.csv"), DataFrame; delim=';');
	scale!(df, [:kcal_per_g, :perc_fat, :perc_lactose])
end;

# ╔═╡ 1c3fc474-fe76-11ea-327f-7beea4fb1872
stan6_4 = "
data{
  int <lower=1> N;              // Sample size
  vector[N] K;
  vector[N] L;
}
parameters{
  real a;
  real bL;
  real<lower=0> sigma;
}
model{
  vector[N] mu;
  sigma ~ exponential( 1 );
  a ~ normal( 0 , 0.2 );
  bL ~ normal( 0 , 0.5 );
  mu = a + bL * L;
  K ~ normal( mu , sigma );
}
";

# ╔═╡ 1c48ffe4-fe76-11ea-0883-f3715e7beda4
md"##### Define the SampleModel, etc."

# ╔═╡ 1c4ba9c4-fe76-11ea-0aff-ef249af5bdb5
begin
	data = (N = size(df, 1), L = df.perc_lactose_s, K = df.kcal_per_g_s)
	init = (a = 0.0, bL = 1.0, sigma = 1.0)
	q6_4s, m6_4s, o6_4s = stan_quap("m6.3", stan6_4; data, init)
	if !isnothing(m6_4s)
		post6_4s_df = read_samples(m6_4s, :dataframe)
		PRECIS(post6_4s_df)
	end
end

# ╔═╡ 1c57c894-fe76-11ea-0b4f-152cfd868c3f
!isnothing(m6_4s) && (part6_4s = Particles(post6_4s_df))

# ╔═╡ 1c586718-fe76-11ea-10ec-892fadc8778c
if !isnothing(q6_4s)
	quap6_4s_df = sample(q6_4s)
	PRECIS(quap6_4s_df)
end

# ╔═╡ 1c659f3c-fe76-11ea-0e95-c1c2cadaaebc
!isnothing(m6_4s) && hpdi(part6_4s.bL.particles, alpha=0.11)

# ╔═╡ 1c6741e8-fe76-11ea-34f3-c54c70b51171
md"## End of clip-06-06-09s.jl"

# ╔═╡ Cell order:
# ╟─d1b31aa0-fe75-11ea-0d96-bdbf19c2c1c7
# ╠═1c31a6c8-fe76-11ea-1122-0580af80a98b
# ╠═1c31e480-fe76-11ea-356d-6b2d21e37ace
# ╠═1c325ef6-fe76-11ea-1967-45317ad16daa
# ╠═1c3fc474-fe76-11ea-327f-7beea4fb1872
# ╟─1c48ffe4-fe76-11ea-0883-f3715e7beda4
# ╠═1c4ba9c4-fe76-11ea-0aff-ef249af5bdb5
# ╠═1c57c894-fe76-11ea-0b4f-152cfd868c3f
# ╠═1c586718-fe76-11ea-10ec-892fadc8778c
# ╠═1c659f3c-fe76-11ea-0e95-c1c2cadaaebc
# ╟─1c6741e8-fe76-11ea-34f3-c54c70b51171


text:
#= Julia syntax highlighting test.

Modified from https://github.com/JuliaEditorSupport/julia-syntax-test-cases

This file derives from https://gist.github.com/Wilfred/f1aca44c61ed6e1df603
whose author is [@Wilfred](https://github.com/Wilfred). @Wilfred has put it in
the public domain:
  https://gist.github.com/Wilfred/f1aca44c61ed6e1df603#gistcomment-2948526

Changes from the original are governed by the license of the repository in
which this file is found.

This file is designed to test various corner cases of Julia
syntax highlighting.
=#

baremodule Mod1
    # Nothing here
end

module Mod2
    # Here neither
end

## Simple function definitions.
# Expected: `function` should be highlighted, as should `foo_bar!`.
function foo_bar!(x,y)
    x + y + 1
end

# Expected: `foo_bar!` should be highlighted.
foo_bar!(x,y) = x + y

# Expected: `foo` should be highlighted.
Baz.foo(x) = 1

# Expected: `foo` should be highlighted.
foo(x::(Int,)) = 1

# Expected: `foo` should be highlighted.
foo(x, y=length(x))

## Function definitions in namespaces.
# Expected: `bar` should be highlighted.
function Foo.bar(x, y)
    x + 1
end

## Function definitions with type variables.
# Expected: `elsize` should be highlighted.
elsize(::AbstractArray{T}) where {T} = sizeof(T)

function elsize(::AbstractArray{T}) where T
    sizeof(T)
end

## Nested brackets in function definitions.
# Expected: `cell` should be highlighted.
cell(dims::(Integer...)) = Array(Any, convert((Int...), dims))

# TODO: find an example with a nested type expression.

## Macro usage
# Expected: `@hello_world!` should be highlighted.
@hello_world! foo

# Expected: highlight `myfun`
@inline myfun() = println("myfun")

## Builtin functions.
# Expected: `throw`, `error` and `super` should not be highlighted. There are
# too many built-in functions for this to be useful.
# https://github.com/JuliaLang/julia/commit/134867c69096fcb52afa5d5a7433892b5127e981
# https://github.com/JuliaLang/julia/pull/7963#issuecomment-52586261
throw(foo)
error("foo", bar, "baz")
super(Int)

## Strings
# Expected: highlight the string.
x = "foo \"bar\" baz"

# Expected: highlight the whole string.
x = "foo
bar"

# Expected: highlight the whole triple-quoted string.
x = """hello "world" foobar"""
y = """foo\\"""
z = """bar\""""
w = """"bar"""

# Expected: highlight `$user`
x = "hello $user"

# Expected: don't highlight `$user`
x = "hello \$user"

# Expected: highlight `$val`
x = """(a="$val")"""

# Expected: treat r as part of the string, so `r"a"` is highlighted.
x = r"0.1"

# Expected: treat ismx as part of the string, so `r"a"ismx` is highlighted.
x = r"a"ismx

# Expected: highlight `r"""a "b" c"""`
x = r"""a "b" c"""

# Expected: treat v as part of the string, so `v"0.1"` is highlighted.
x = v"0.1"

# Expected: treat b as part of the string, so `b"a"` is highlighted.
x = b"a"

# Bonus points:
# Expected: highlight the interpolation brackets `$(` and `)`
x = "hello $(user * user)"

# Bonus points:
# Expected: highlight regexp metacharacters `[` and `]`
x = r"[abc]"

# Bonus points:
# Expected: highlight escape sequences `\xff` and `\u2200`
x = b"DATA\xff\u2200"

# Bonus points:
# Expected: don't highlight `$user`
x = raw"hello $user"

## Characters
# Expected: highlight the character.
x = 'a'
y = '\u0'
z = '\U10ffff'
w = '\x41'
a = ' '
b = '"'
c = '''
d = '\''
e = '\\'

# Expected: don't highlight, as ' is an operator here, not a character delimiter.
a = b' + c'
A'''

# Bonus points:
# Expected: don't highlight the character
x = 'way too long so not a character'
x = ''

## Comments
# Expected: highlight `# foo`
# foo

# Expected: highlight `#= foo\n bar =#`
#= foo
bar =#

# Expected: highlight `#= #= =# =#` (comments can nest).
#= #= =# =#

# Expected: highlight `'` as adjoint operator
A#==#'
(A)#==#'
A[1]#==#'

## Type declarations

# Expected highlight `Foo` and `Bar`
mutable struct Foo
    x::Bar
end

# Expected highlight `Foo` and `Bar`
struct Foo
    x::Bar
end

# Expected: highlight `Foo` and `Bar`
abstract type Foo <: Bar end

# Expected: don't highlight x or y
x <: y

## Type annotations

# Expected: highlight `FooBar`
f(x::FooBar) = x

# Expected: highlight `Int8`
function foo()
    local x::Int8 = 5
    x
end

# Expected: highlight `Point` and `Real` as types
function norm(p::Point{<:Real})
    sqrt(p.x^2 + p.y^2)
end

# Expected: highlight `g` as function and `Int8` as type
function g(x, y)::Int8
   return x * y
end

# Expected: highlight `T` and `Number`
same_type_numeric(x::T, y::T) where {T <: Number} = true
same_type_numeric(x::T, y::T) where T = false

## Parametric type declaration

# Expected: highlight `Pointy` and `T`
abstract type Pointy{T} end

# Expected: highlight `Point`, `Pointy` and `T`
struct Point{T} <: Pointy{T}
    x::T
    y::T
end

## Variable declarations

# Expected: highlight `x` and `y`
global x = "foo, bar = 2", y = 3

# Expected: highlight `x` and `y`
global x = foo(a, b), y = 3

# Expected: highlight `y`
const y = "hello world"

# Expected: highlight `x` and `y`
function foo()
    local x = f(1, 2), y = f(3, 4)
    x + y
end

# Expected: highlight `x` and `y`
let x = f(1, 2), y = f(3, 4)
    x + y
end

## Colons and end

# Expected: highlight `:foo`, `:end` and `:function`
:foo
x = :foo
y = :function
z = :end

# Expected: highlight index `[end]` differently to block delimiter `end`
if foo[end]
end

# Expected: highlight as index `end`
foo[bar:end]

# Expected: highlight as index `begin`
foo[begin:4]

# Expected: don't highlight `:123`
x = :123

# Expected: don't highlight `:baz`
foo[bar:baz]

# Expected: highlight `:baz`
foo[:baz]

# Expected: highlight both `:baz`
foo(:baz,:baz)

# Note that `: foo` is currently a valid quoted symbol, this will hopefully
# change in 0.4: https://github.com/JuliaLang/julia/issues/5997

# Expected: highlight `:foo`
[1 :foo]

# Expected: highlight `:end`
[1 :end]

# Expected: highlight `:two`
@eval :one+:two

# Expected: don't highlight `:end` but `end` as index
[(1+1):end]

# Expected: don't highlight `:end` but `end` as index
[a[1]:end]

# Expected: don't highlight `:foo`
for x=1:foo
    print(x)
end

## Range detection

# Bonus points:
# Expected: don't highlight `:s2`
push!(a, s1 :s2)

# Bonus points:
# Expected: don't highlight `:end`
a[begin :end]

## Expression evaluation

# Expected: highlight `:` as operator
a = :(x = 2)

# Expected: highlight `:call` and `:b` as symbols
# Debatable: highlight `:+` as operator
ex = Expr(:call, :+, a, :b)

## Number highlighting

# Expected: highlight all these as numbers
x = 123
x = 1.1
x = .5
x = 0x123abcdef
x = 0o7
x = 0b1011
x = 2.5e-4
x = 2.5E-4
x = 1e+00
x = 2.5f-4
x = 0x.4p-1
x = 1_000

# Expected: highlight these as numbers or built-ins
x = Inf
x = NaN

# Expected: highlight `123`, not the letter
y = 123x
y = 123e

# Expected: highlight `1e+1` and `1e-1`
1e+1+1e-1

# Expected: highlight `1.` and `.1`
1. +.1
# Note that `1.+1` is currently ambiguous and gives an error

# Bonus points:
# Expected: don't highlight `..`
x = 1..3

# Bonus points:
# Debatable: highlight the first two digits, not the following digits
# or show an error
x = 0o1291
x = 0b1091

# Debatable: highlight `π` as a number or built-in
# (note that `πx` is a single symbol, not `π * x`)
x = π

## List comprehension
# Expected: highlight `for` and `if` without the `end` keyword
[i for i in 1:5 if i%2==0]

## Broadcasting
# Expected: highlight `.+` as operator
a.+1

## Command
# Expected: highlight "`echo 1`"
c = `echo 1`

# Expected: highlight "```echo `hello 1` ```"
c = ```echo `hello 1` ```

# Expected: highlight "raw`echo $1`"
c = raw`echo $1`

## Non-standard identifiers
# Bonus points:
# Expected: highlight `var"##"` as a function
function var"##"(x)
    println(x)
end

# Bonus points:
# Expected: highlight `var"%%"` as a function
var"%%"(x) = println(x)

# Bonus points:
# Expected: highlight `$var` as string and `##""` as comment
"$var"##""

# Bonus points:
# Expected: highlight `$(var")(")` as string interpolation
"$(var")(")"

# Bonus points:
# Expected: highlight `'` as adjoint operator
var"##mat"'

## Code folding: for and if in list comprehension
# Expected: fold between function and last end
function test(x)
    a = (if x; 0 else 1 end)
    println(a)
end


text:
using Test
using Profile: Allocs

@testset "alloc profiler doesn't segfault" begin
    res = Allocs.@profile sample_rate=1.0 begin
        # test the allocations during compilation
        using Base64
    end
    profile = Allocs.fetch()

    @test length(profile.allocs) > 0
    first_alloc = profile.allocs[1]
    @test first_alloc.size > 0
    @test length(first_alloc.stacktrace) > 0
    @test length(string(first_alloc.type)) > 0
end

@testset "alloc profiler works when there are multiple tasks on multiple threads" begin
    NUM_TASKS = 1000

    # This test is only really meaningful if we're running on
    # multiple threads, but this isn't true on the windows tests,
    # causing them to fail. So, commenting this assertion out.
    # @test Threads.nthreads() > 1

    function do_work()
        ch = Channel{Vector{Float64}}(Inf)
        @sync for i in 1:NUM_TASKS
            Threads.@spawn begin
                # generate garbage
                put!(ch, zeros(100))
            end
        end
        close(ch)
    end

    # call once to make sure it's compiled
    precompile(do_work, ())
    do_work()

    res = Allocs.@profile sample_rate=1 begin
        do_work()
    end
    profile = Allocs.fetch()

    # expecting at least 2 allocations per task:
    # 1. the task
    # 2. the vector
    @test length(profile.allocs) >= 2*NUM_TASKS
    first_alloc = profile.allocs[1]
    @test first_alloc.size > 0
    @test length(first_alloc.stacktrace) > 0
    @test length(string(first_alloc.type)) > 0

    @testset for type in (Task, Vector{Float64},)
        @test length(filter(a->a.type <: type, profile.allocs)) >= NUM_TASKS
    end

    # TODO: it would be nice to assert that these tasks
    # were actually scheduled onto multiple threads,
    # and we see allocs from all threads in the profile
end

@testset "alloc profiler start stop fetch clear" begin
    function do_work()
        # Compiling allocates a lot
        for f in (gensym() for _ in 1:10)
            @eval begin
                $f() = 10
                $f()
            end
        end
    end

    Allocs.@profile sample_rate=1 do_work()
    @test length(Allocs.fetch().allocs) > 10

    Allocs.clear()
    @test length(Allocs.fetch().allocs) == 0
    Allocs.clear()
    @test length(Allocs.fetch().allocs) == 0

    Allocs.@profile sample_rate=1 do_work()
    curr_allocs = length(Allocs.fetch().allocs)
    @test curr_allocs > 10

    # Do _more_ work, adding into the same profile
    Allocs.@profile sample_rate=1 do_work()
    @test length(Allocs.fetch().allocs) > curr_allocs

    Allocs.clear()
    @test length(Allocs.fetch().allocs) == 0

    # Clear without fetching

    Allocs.@profile sample_rate=1 do_work()
    Allocs.clear()
    @test length(Allocs.fetch().allocs) == 0

    # And things still work like normal afterwards

    Allocs.@profile sample_rate=1 do_work()
    Allocs.@profile sample_rate=1 do_work()
    Allocs.@profile sample_rate=1 do_work()
    @test length(Allocs.fetch().allocs) > 10

    Allocs.@profile sample_rate=1 do_work()
    Allocs.@profile sample_rate=1 do_work()
    @test length(Allocs.fetch().allocs) > 10

    Allocs.clear()
end

@testset "alloc profiler catches strings" begin
    Allocs.@profile sample_rate=1 "$(rand())"

    prof = Allocs.fetch()
    Allocs.clear()

    @test length(prof.allocs) >= 1
    @test length([a for a in prof.allocs if a.type == String]) >= 1
end

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia


text:
@testset "Traits" begin
  @testset "Domain" begin
    # dummy type implementing the Domain trait
    struct DummyDomain{Dim,T} <: Domain{Dim,T}
      origin::Point{Dim,T}
    end
    function Base.getindex(domain::DummyDomain{Dim,T}, ind::Int) where {Dim,T}
      c = domain.origin + Vec(ntuple(i->T(ind), Dim))
      r = one(T)
      Ball(c, r)
    end
    Meshes.nelements(d::DummyDomain) = 3

    # basic properties
    dom = DummyDomain(P2(0,0))
    @test embeddim(dom) == 2
    @test coordtype(dom) == T
    @test eltype(dom) <: Ball{2,T}

    # coordinates of centroids
    dom = DummyDomain(P2(1,1))
    @test coordinates(dom, 1:3) == T[2 3 4; 2 3 4]

    dom = DummyDomain(P2(0,0))
    @test sprint(show, dom) == "3 DummyDomain{2,$T}"
    if T == Float32
      @test sprint(show, MIME"text/plain"(), dom) == "3 DummyDomain{2,Float32}\n  └─Ball{2,Float32}(Point(1.0f0, 1.0f0), 1.0))\n  └─Ball{2,Float32}(Point(2.0f0, 2.0f0), 1.0))\n  └─Ball{2,Float32}(Point(3.0f0, 3.0f0), 1.0))"
    elseif T == Float64
      @test sprint(show, MIME"text/plain"(), dom) == "3 DummyDomain{2,Float64}\n  └─Ball{2,Float64}(Point(1.0, 1.0), 1.0))\n  └─Ball{2,Float64}(Point(2.0, 2.0), 1.0))\n  └─Ball{2,Float64}(Point(3.0, 3.0), 1.0))"
    end

    if visualtests
      dom = DummyDomain(P2(0,0))
      @test_reference "data/domain-$T.png" plot(dom)
      @test_reference "data/domain-data-$T.png" plot(dom,1:3)
    end
  end

  @testset "Data" begin
    # dummy type implementing the Data trait
    struct DummyData{𝒟,𝒯} <: Data
      domain::𝒟
      table::𝒯
    end
    Meshes.domain(data::DummyData) = data.domain
    Meshes.values(data::DummyData) = data.table
    Meshes.constructor(::Type{D}) where {D<:DummyData} = DummyData

    # equality of data sets
    data₁ = DummyData(CartesianGrid{T}(2,2), (a=[1,2,3,4], b=[5,6,7,8]))
    data₂ = DummyData(CartesianGrid{T}(2,2), (a=[1,2,3,4], b=[5,6,7,8]))
    data₃ = DummyData(PointSet(rand(P2,4)), (a=[1,2,3,4], b=[5,6,7,8]))
    @test data₁ == data₂
    @test data₁ != data₃
    @test data₂ != data₃

    # Tables interface
    dom = CartesianGrid{T}(2,2)
    dat = DummyData(dom, (a=[1,2,3,4], b=[5,6,7,8]))
    @test Tables.istable(dat)
    @test Tables.rowaccess(dat)
    s = Tables.schema(dat)
    @test s.names == (:a,:b,:geometry)
    @test s.types == (Int, Int, Quadrangle{2,T,Vector{P2}})
    @test collect(Tables.rows(dat)) == [
      (1, 5, dom[1]),
      (2, 6, dom[2]),
      (3, 7, dom[3]),
      (4, 8, dom[4])
    ]

    # variables interface
    data = DummyData(CartesianGrid{T}(2,2), (a=[1,2,3,4], b=[5,6,7,8]))
    @test data[:a] == [1,2,3,4]
    @test data[:b] == [5,6,7,8]
    @test_throws ErrorException data[:c]
    @test_throws ErrorException data[:geometry]

    # utility functions
    data = DummyData(PointSet(rand(P2,4)), (a=[1,2,3,4], b=[5,6,7,8]))
    @test asarray(data, :a) == [1,2,3,4]
    @test asarray(data, :b) == [5,6,7,8]
    data = DummyData(CartesianGrid{T}(2,2), (a=[1,2,3,4], b=[5,6,7,8]))
    @test asarray(data, :a) == [1 3; 2 4]
    @test asarray(data, :b) == [5 7; 6 8]

    data = DummyData(CartesianGrid{T}(2,2), (a=[1,2,3,4], b=[5,6,7,8]))
    @test sprint(show, data) == "4 DummyData{2,$T}"
    if T == Float32
      @test sprint(show, MIME"text/plain"(), data) == "4 DummyData{2,Float32}\n  variables\n    └─a (Int64)\n    └─b (Int64)\n  domain: 2×2 CartesianGrid{2,Float32}"
    elseif T == Float64
      @test sprint(show, MIME"text/plain"(), data) == "4 DummyData{2,Float64}\n  variables\n    └─a (Int64)\n    └─b (Int64)\n  domain: 2×2 CartesianGrid{2,Float64}"
    end

    if visualtests
      data = DummyData(CartesianGrid{T}(2,2), (a=[1,2,3,4], b=[5,6,7,8]))
      @test_reference "data/data-$T.png" plot(data)
    end
  end
end


text:
# # Pyramid example
#
# This is the Pyramid example in the book "Evaluate Derivatives", Sec. 3.5.

using NiLang, NiLang.AD

@i function pyramid!(y!, v!, x::AbstractVector{T}) where T
    @safe @assert size(v!,2) == size(v!,1) == length(x)
    @invcheckoff @inbounds for j=1:length(x)
        v![1,j] += x[j]
    end
    @invcheckoff @inbounds for i=1:size(v!,1)-1
        for j=1:size(v!,2)-i
            @routine begin
                @zeros T c s
                c += cos(v![i,j+1])
                s += sin(v![i,j])
            end
            v![i+1,j] += c * s
            ~@routine
        end
    end
    y! += v![end,1]
end

x = randn(20)
pyramid!(0.0, zeros(20, 20), x)

# Let' compare our implementation with those in the book
# ![](pyramid-benchmark.png)
using BenchmarkTools
@benchmark gradient(Val(1), pyramid!, (0.0, zeros(20, 20), $x))


text:
using Knet
include("readdata.jl")

# Try a convolutional model, similar to CharCNN in http://arxiv.org/pdf/1508.06615.pdf
# The input is a sequence of characters of length T.
# We map this to an embedding matrix DxT (alternative is to use one-hot with D=4).
# We apply a 1-D convolution bank (width and number of filters to be decided).
# In the simplest case with one conv layer, we apply max pooling of window T.
# We apply an affine transform and a softmax to the output.
# B shows the minibatch size, and everything is in 4D for CUDNN ops.

# x is a (1,T,D,N) matrix of inputs where T is the sequence length, D is the dimension of each char embedding, N is the minibatch size
# we convolve with a (1,W,D,C) filter bank with W/2 padding to avoid dimension change
# this gives a (1,T,C,N) output which gets a bias added and f applied, neither of which effects the dimension
# we apply max pooling with pwindow=T which gives an output of (1,1,C,N)
# we have a fully connected layer ending with a softmax at the end
# cbfp73 has the following default options each of which can be modified in model01:
# cwindow=3, cpadding=div(cwindow,2), cstride=1, cmode=Knet.CUDNN_CONVOLUTION,
# pwindow=2, ppadding=0, pstride=pwindow, pmode=Knet.CUDNN_POOLING_MAX,
# cinit=Xavier(), binit=Constant(0), f=:relu

@knet function model01(x; T=51, D=4, C=512, W=3, o...)
    y = cbfp73(x; o..., pwindow=T, out=C)
    return wbf73(y; o..., out=2, f=:soft)
end

function train01(; N=2^13, lr=0.001, adam=true, nbatch=128, o...)
    global f = compile(:model01; o...)
    setp(f, lr=lr, adam=adam)
    global d = GulayseData(nbatch)
    sloss = zloss = -1; n = nextn = 1
    for (x,ygold) in d
        ypred = forw(f, x)
        sl = softloss(ypred,ygold); sloss = (n==1 ? sl : 0.99 * sloss + 0.01 * sl)
        zl = zeroone(ypred,ygold);  zloss = (n==1 ? zl : 0.99 * zloss + 0.01 * zl)
        n==nextn && (println((n,sloss,1-zloss)); nextn*=2)
        back(f, ygold, softloss)
        update!(f)
        (n += 1) > N && break
    end
end

train01()

max_stars_repo_name: JuliaTagBot/BindingPrediction

max_issues_repo_name: JuliaTagBot/BindingPrediction

max_forks_repo_name: JuliaTagBot/BindingPrediction


text:
using Base.Test
include("reduce.jl")

dev = CuDevice(0)
if capability(dev) < v"3.0"
    warn("this example requires a newer GPU")
    exit(0)
end

len = 10^7
input = ones(Int32, len)

# CPU
cpu_val = reduce(+, input)

# CUDAnative
let
    ctx = CuContext(dev)
    gpu_input = CuArray(input)
    gpu_output = similar(gpu_input)
    gpu_reduce(+, gpu_input, gpu_output)
    gpu_val = Array(gpu_output)[1]
    destroy(ctx)
    @assert cpu_val == gpu_val
end


text:
@doc """
    S = read_hdf5(filestr, s::TimeSpec, t::TimeSpec, [, keywords])
    read_hdf5!(S, filestr, s::TimeSpec, t::TimeSpec, [, keywords])

Read data in seismic HDF5 format from files matching pattern `filestr`.

`s`, `t` are required arguments but can be any Type ∈ (DateTime, Real, String);
type `?TimeSpec` for more information about how these are interpreted.

|KW     | Type      | Default   | Meaning                                     |
|:---   |:---       |:---       |:---                                         |
| id    | String    | "*"       | id pattern, formated nn.sss.ll.ccc          |
|       |           |           |  (net.sta.loc.cha); FDSN wildcards [^1]     |
| msr   | Bool      | true      | read full (MultiStageResp) instrument resp? |
| v     | Integer   | 0         | verbosity                                   |

[^1] A question mark ('?') is a wildcard for a single character; an asterisk ('*') is a wildcard for zero or more characters

See also: `TimeSpec`, `parsetimewin`, `read_data`
""" read_hdf5!
function read_hdf5!(S::GphysData, fpat::String, s::TimeSpec, t::TimeSpec;
  fmt ::String                = "asdf",                 # data format
  id  ::Union{String, Regex}  = "*",                    # id string
  msr ::Bool                  = true,                   # read multistage response?
  v   ::Integer               = KW.v                    # verbosity
  )

  N = S.n
  filestr = abspath(fpat)
  one_file = safe_isfile(filestr)

  if fmt == "asdf"
    if one_file
      append!(S, read_asdf(filestr, id, s, t, msr, v))
    else
      files = ls(filestr)
      for fname in files
        append!(S, read_asdf(fname, id, s, t, msr, v))
      end
    end

  else
    error("Unknown file format (possibly NYI)!")
  end

  new_chan_src = view(S.src, N+1:S.n)
  fill!(new_chan_src, filestr)
  note!(S, N+1:S.n, string( " ¦ +source ¦ read_hdf5!(S, ",
                            "\"", fmt, "\", ",
                            "\"", s, "\", ",
                            "\"", t, "\", ",
                            "fmt=\"", fmt, "\", ",
                            "id=\"", id, "\", ",
                            "msr=", msr, ", ",
                            "v=", KW.v, ")")
        )

  return nothing
end

@doc (@doc read_hdf5!)
function read_hdf5(filestr::String, s::TimeSpec, t::TimeSpec;
  fmt ::String                = "asdf",                 # data format
  id  ::Union{String, Regex}  = "*",                    # id string
  msr ::Bool                  = true,                   # read multistage response?
  v   ::Integer               = KW.v                    # verbosity
  )

  S = SeisData()
  read_hdf5!(S, filestr, s, t,
    fmt     = fmt,
    id      = id,
    msr     = msr,
    v       = v
    )
  return S
end


text:
import Base: -, +, *, /, ^, abs, abs2, sqrt, exp, exp2, expm1, exp10, log, log2,
    log1p, log10, sin, sinpi, cos, cospi, tan, asin, acos, atan, sinh, cosh,
    asinh, acosh, tanh, atanh, inv, floor, ceil, min, max, round, trunc, eps
export -, +, *, /, ^, abs, abs2, sqrt, exp, exp2, expm1, exp10, log, log2,
    log1p, log10, sin, sinpi, cos, cospi, tan, asin, acos, atan, sinh, cosh,
    asinh, acosh, tanh, atanh, inv, floor, ceil, min, max, round, trunc, eps

for f in (:-, :abs, :abs2, :sqrt, :exp, :exp2, :expm1, :exp10, :log, :log2,
          :log1p, :log10, :sin, :sinpi, :cos, :cospi, :tan, :asin, :acos, :atan,
          :sinh, :cosh, :asinh, :acosh, :tanh, :atanh, :inv, :floor, :ceil,
          :round, :trunc)
    @eval $f(a::NumberInterval) = NumberInterval($f(Interval(a)))
end

for f in (:+, :-, :*, :/, :^, :atan, :min, :max)
    @eval $f(a::NumberInterval, b::NumberInterval) =
        NumberInterval($f(Interval(a), Interval(b)))
end

# specialization necessary to avoid ambiguity
for t in (:Integer, :AbstractFloat, :Rational)
    @eval ^(a::NumberInterval, n::$t) = NumberInterval(Interval(a)^n)
end

round(a::NumberInterval, mode) = NumberInterval(round(Interval(a), mode))

eps(::Type{NumberInterval{T}}) where T = eps(T)
function eps(a::NumberInterval)
    _a = abs(a)
    return NumberInterval(eps(_a.lo), eps(_a.hi))
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    ColTable()

The transform that applies the function `Tables.columntable` to to the input table.
"""
struct ColTable <: Stateless end

apply(::ColTable, table) = Tables.columntable(table), table

revert(::ColTable, newtable, cache) = cache


text:
using Histograms
using Base.Test

h1 = ErrorHistogram(linspace(0, 1, 11))
@test ndims(h1) == 1
@test nbins(h1) == (10, )

h2 = ErrorHistogram(linspace(0, 1, 11), linspace(0, 1, 3))
@test ndims(h2) == 2
@test nbins(h2) == (10, 2)

push!(h1, 0.0)
push!(h1, 0.0)
push!(h1, 0.2, 0.5)
push!(h1, 0.2, 0.5)
push!(h1, 0.2, 0.5)

@test all(contents(h1) .== Float64[2.0, 0.0, 1.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
@test_approx_eq sum(errors(h1) .- Float64[sqrt(2.0), 0.0, sqrt(0.75), 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) 0.0

push!(h2, (0.0, 0.8), 0.5)
@test contents(h2)[1,2] == 0.5

htot = h2 + h2

@test_approx_eq integral(htot)[1] 2*integral(h2)[1]


text:
using MPI
using Test
using LinearAlgebra
using Random
using StaticArrays
using ClimateMachine
using ClimateMachine.SystemSolvers
using ClimateMachine.MPIStateArrays
using CUDA
using Random
using KernelAbstractions

import ClimateMachine.MPIStateArrays: array_device

ClimateMachine.init(; fix_rng_seed = true)

@kernel function multiply_by_A!(x, A, y, n1, n2)
    I = @index(Global)
    for i in 1:n1
        tmp = zero(eltype(x))
        for j in 1:n2
            tmp += A[i, j, I] * y[j, I]
        end
        x[i, I] = tmp
    end
end

let
    if CUDA.has_cuda_gpu()
        Arrays = [Array, CuArray]
    else
        Arrays = [Array]
    end

    for ArrayType in Arrays
        for T in [Float32, Float64]
            ϵ = eps(T)

            @testset "($ArrayType, $T) Basic Test" begin
                Random.seed!(42)

                # Test 1: Basic Functionality
                n = 100   # size of local (batch) matrix
                ni = 100  # batch size

                b = ArrayType(randn(n, ni))  # rhs
                x = ArrayType(randn(n, ni))  # initial guess
                x_ref = similar(x)

                A = ArrayType(randn((n, n, ni)) ./ sqrt(n))
                for i in 1:n
                    A[i, i, :] .+= 10i
                end

                ss = size(b)[1]
                bgmres = BatchedGeneralizedMinimalResidual(
                    b,
                    n,
                    ni,
                    M = ss,
                    atol = ϵ,
                    rtol = ϵ,
                )

                # Define the linear operator
                function closure_linear_operator_multi!(A, n1, n2, n3)
                    function linear_operator!(x, y)
                        device = array_device(x)
                        if isa(device, CPU)
                            groupsize = Threads.nthreads()
                        else # isa(device, CUDADevice)
                            groupsize = 256
                        end
                        event = Event(device)
                        event = multiply_by_A!(device, groupsize)(
                            x,
                            A,
                            y,
                            n1,
                            n2,
                            ndrange = n3,
                            dependencies = (event,),
                        )
                        wait(device, event)
                        nothing
                    end
                end
                linear_operator! = closure_linear_operator_multi!(A, size(A)...)

                # Now solve
                linearsolve!(
                    linear_operator!,
                    nothing,
                    bgmres,
                    x,
                    b;
                    max_iters = Inf,
                )

                # reference solution
                for i in 1:ni
                    x_ref[:, i] = A[:, :, i] \ b[:, i]
                end

                @test norm(x - x_ref) < 1000ϵ
            end

            ###
            # Test 2: MPI State Array
            ###
            @testset "($ArrayType, $T) MPIStateArray Test" begin
                Random.seed!(43)

                n1 = 8
                n2 = 3
                n3 = 10
                mpicomm = MPI.COMM_WORLD
                mpi_b = MPIStateArray{T}(mpicomm, ArrayType, n1, n2, n3)
                mpi_x = MPIStateArray{T}(mpicomm, ArrayType, n1, n2, n3)
                mpi_A = ArrayType(randn(n1 * n2, n1 * n2, n3))

                mpi_b.data[:] .= ArrayType(randn(n1 * n2 * n3))
                mpi_x.data[:] .= ArrayType(randn(n1 * n2 * n3))

                bgmres = BatchedGeneralizedMinimalResidual(
                    mpi_b,
                    n1 * n2,
                    n3,
                    M = n1 * n2,
                    atol = ϵ,
                    rtol = ϵ,
                )

                # Now define the linear operator
                function closure_linear_operator_mpi!(A, n1, n2, n3)
                    function linear_operator!(x, y)
                        alias_x = reshape(x.data, (n1, n3))
                        alias_y = reshape(y.data, (n1, n3))
                        device = array_device(x)
                        if isa(device, CPU)
                            groupsize = Threads.nthreads()
                        else # isa(device, CUDADevice)
                            groupsize = 256
                        end
                        event = Event(device)
                        event = multiply_by_A!(device, groupsize)(
                            alias_x,
                            A,
                            alias_y,
                            n1,
                            n2,
                            ndrange = n3,
                            dependencies = (event,),
                        )
                        wait(device, event)
                        nothing
                    end
                end
                linear_operator! =
                    closure_linear_operator_mpi!(mpi_A, size(mpi_A)...)

                # Now solve
                linearsolve!(
                    linear_operator!,
                    nothing,
                    bgmres,
                    mpi_x,
                    mpi_b;
                    max_iters = Inf,
                )

                # check all solutions
                norms = -zeros(n3)
                for cidx in 1:n3
                    sol =
                        Array(mpi_A[:, :, cidx]) \
                        Array(mpi_b.data[:, :, cidx])[:]
                    norms[cidx] = norm(sol - Array(mpi_x.data[:, :, cidx])[:])
                end
                @test maximum(norms) < 8000ϵ
            end

            ###
            # Test 3: Columnwise test
            ###
            @testset "(Array, $T) Columnwise Test" begin

                Random.seed!(2424)
                function closure_linear_operator_columwise!(A, tup)
                    function linear_operator!(y, x)
                        alias_x = reshape(x, tup)
                        alias_y = reshape(y, tup)
                        for i6 in 1:tup[6]
                            for i4 in 1:tup[4]
                                for i2 in 1:tup[2]
                                    for i1 in 1:tup[1]
                                        tmp = alias_x[i1, i2, :, i4, :, i6][:]
                                        tmp2 = A[i1, i2, i4, i6] * tmp
                                        alias_y[i1, i2, :, i4, :, i6] .=
                                            reshape(tmp2, (tup[3], tup[5]))
                                    end
                                end
                            end
                        end
                    end
                end

                tup = (3, 4, 7, 6, 5, 2)
                B = [
                    randn(tup[3] * tup[5], tup[3] * tup[5])
                    for
                    i1 in 1:tup[1],
                    i2 in 1:tup[2], i4 in 1:tup[4], i6 in 1:tup[6]
                ]
                columnwise_A = [
                    B[i1, i2, i4, i6] + 10I
                    for
                    i1 in 1:tup[1],
                    i2 in 1:tup[2], i4 in 1:tup[4], i6 in 1:tup[6]
                ]
                # taking the inverse of A isn't great, but it is convenient
                columnwise_inv_A = [
                    inv(columnwise_A[i1, i2, i4, i6])
                    for
                    i1 in 1:tup[1],
                    i2 in 1:tup[2], i4 in 1:tup[4], i6 in 1:tup[6]
                ]
                columnwise_linear_operator! =
                    closure_linear_operator_columwise!(columnwise_A, tup)
                columnwise_inverse_linear_operator! =
                    closure_linear_operator_columwise!(columnwise_inv_A, tup)

                mpi_tup = (tup[1] * tup[2] * tup[3], tup[4], tup[5] * tup[6])
                b = randn(mpi_tup)
                x = copy(b)
                columnwise_inverse_linear_operator!(x, b)
                x +=
                    randn((tup[1] * tup[2] * tup[3], tup[4], tup[5] * tup[6])) *
                    0.1

                reshape_tuple_f = tup
                permute_tuple = (5, 3, 1, 4, 2, 6)

                bgmres = BatchedGeneralizedMinimalResidual(
                    b,
                    tup[3] * tup[5],
                    tup[1] * tup[2] * tup[4] * tup[6];
                    M = tup[3] * tup[5],
                    forward_reshape = tup,
                    forward_permute = permute_tuple,
                    atol = 10ϵ,
                    rtol = 10ϵ,
                )

                x_exact = copy(x)
                linearsolve!(
                    columnwise_linear_operator!,
                    nothing,
                    bgmres,
                    x,
                    b,
                    max_iters = tup[3] * tup[5],
                )

                columnwise_inverse_linear_operator!(x_exact, b)
                @test norm(x - x_exact) / norm(x_exact) < 1000ϵ
                columnwise_linear_operator!(x_exact, x)
                @test norm(x_exact - b) / norm(b) < 1000ϵ
            end
        end
    end
end


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

mutable struct IoK8sApiCoreV1TCPSocketAction <: SwaggerModel
    host::Any # spec type: Union{ Nothing, String } # spec name: host
    port::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: port

    function IoK8sApiCoreV1TCPSocketAction(;host=nothing, port=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1TCPSocketAction, Symbol("host"), host)
        setfield!(o, Symbol("host"), host)
        validate_property(IoK8sApiCoreV1TCPSocketAction, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        o
    end
end # type IoK8sApiCoreV1TCPSocketAction

const _property_map_IoK8sApiCoreV1TCPSocketAction = Dict{Symbol,Symbol}(Symbol("host")=>Symbol("host"), Symbol("port")=>Symbol("port"))
const _property_types_IoK8sApiCoreV1TCPSocketAction = Dict{Symbol,String}(Symbol("host")=>"String", Symbol("port")=>"IoK8sApimachineryPkgUtilIntstrIntOrString")
Base.propertynames(::Type{ IoK8sApiCoreV1TCPSocketAction }) = collect(keys(_property_map_IoK8sApiCoreV1TCPSocketAction))
Swagger.property_type(::Type{ IoK8sApiCoreV1TCPSocketAction }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1TCPSocketAction[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1TCPSocketAction }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1TCPSocketAction[property_name]

function check_required(o::IoK8sApiCoreV1TCPSocketAction)
    (getproperty(o, Symbol("port")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1TCPSocketAction }, name::Symbol, val)
end