text:
#README
#calc_ttvfast takes a planetary system, plots TTV vs transit time, and returns an array of arrays holding the TTVs for each planet. It is slower than calc_ttvfaster but works for more systems.
#calc_ttvfaster takes a planetary system, plots TTV vs transit time, and returns an array of arrays holding the TTVs for each planet.
#plot_ttv takes an array of TTVs, an array of transit times, and an array of periods for each planet in a system, and plots TTVs vs. transit_times.
#plot_ttv takes an array of TTVs, an array of periods, and an array of first-transit-times (t0) for each planet in a system, and creates an array of transit times before calling plot_ttvfast.
#get_ttv_sum_stats takes an array of TTVs for each planet in a system and returns the scatter, rms and iqr for all the ttvs associated with each planet.

#This code was tested outside of SysSim, to incorporate into SysSim, put the following lines in a "calc_transit_obs" function just after the for loop over each planet.
#    if num_planets(sys)>1
#      sysTTVs = calc_ttvfaster(sys) #,sysTTVs,sys_stats)		#returns an array of arrays of TTVs for each planet calculated.
#      sysTTVs = calc_ttvfast(sys) #,sysTTVs,sys_stats)               	#returns an array of arrays of TTVs for each planet calculated.
#      ttv_sum_stat = get_ttv_sum_stats(sysTTVs)			#returns the scatter, rms and iqr for all the ttvs associated with each planet.
#    end

include("TTVFaster/src/TTVFaster.jl")
include("TTVFast/src/TTVFast.jl")
#include("Stats/src/Stats.jl") required for calculating iqr but the package is currently out of date.
using PyPlot
using StatsBase

# Integration parameters
inflag = 0
t_start=0.0
t_stop=4*365.2425
dt = 0.02
kep_obs_duration = t_stop-t_start
scatter_const = 1.4826;

function calc_ttvfast(sys::Any) #sys::ExoplanetsSysSim.PlanetarySystem{ExoplanetsSysSim.Star}
#calculates transit timing variations for each planet in a system over the course of Kepler's observations.
  nplanets = num_planets(sys)
  p = Array{Cdouble}(undef,2+nplanets*7)		#array holding information for each planet in the system
  p[1] = 4pi/365.2425^2; 				# G
  p[2] = 1.0; 						# Mstar
  num_events = 0;
  for i in 1:nplanets
    period = sys.orbit[i].P
    p[2+7*(i-1)+1] = sys.planet[i].mass 		# Planet Mass
    p[2+7*(i-1)+2] = period 	 			# Period
    p[2+7*(i-1)+3] = sys.orbit[i].ecc			# Eccentricity
    p[2+7*(i-1)+4] = sys.orbit[i].incl*180/pi		# Inclination
    p[2+7*(i-1)+5] = sys.orbit[i].asc_node*180/pi 	# Longitude of Ascending Node
    p[2+7*(i-1)+6] = sys.orbit[i].omega*180/pi		# Argument of Pericenter
    p[2+7*(i-1)+7] = sys.orbit[i].mean_anom*180/pi	# Mean Anomaly
    num_events += ceil(Integer, kep_obs_duration/p[2+7*(i-1)+2] +1); # /* large enough to fit all the transits calculated by the code*/
  end
  ttvfast_input = TTVFast.ttvfast_inputs_type(p, t_start=t_start, t_stop=t_stop, dt=dt)
  incl_rvs = false #if true, Make sure first rv_time is _after_ t_start+dt/2 or else won't get any RV outputs
  if incl_rvs
    num_rvs = 100
    rv_times = collect(range(t_start+0.501*dt, stop=t_stop, length=num_rvs))
    ttvfast_output = TTVFast.ttvfast_outputs_type(num_events ,rv_times)
  else
    ttvfast_output = TTVFast.ttvfast_outputs_type(num_events)
  end
  println(stderr, "# About to call TTVFast")
  @time TTVFast.ttvfast!(ttvfast_input,ttvfast_output)
  sysTTVs = Array{Array{Float64,1}}(undef,nplanets)   		#an array of all the TTVs for each planet
  transit_times = Array{Array{Float64}}(undef,nplanets)		#arrays for each planet holding actual transit time as measured from the beginning of kepler observation.
  ttv_per = zeros(nplanets)					#stores period values for plotting
  for k in 1:nplanets
    sysTTVs[k] = [0.0]
    transit_times[k]=[0.0]
  end
  for j in 1:num_events						#extract transit times
    this_planet = TTVFast.get_event(ttvfast_output,j).planet+1	#determine what planet we are extracting ttvs for
    this_time = TTVFast.get_event(ttvfast_output,j).time	#extract ttv time
    if this_time != 0.0						#events marked 0.0 are extras
      push!(transit_times[this_planet],this_time)
    end
  end
  slope = 0.0
  intercept = 0.0
  for k in 1:nplanets
    popfirst!(transit_times[k])					#remove initializing 0.0 values
    num_times = length(transit_times[k])
    #perform linear regression on transit_times to determine actual period (slope) and t0 (intercept)
    xmean = (num_times-1)/2
    ymean = StatsBase.mean(transit_times[k])
    numer = 0.0
    denom = 0.0
    for m in 1:num_times
      numer += ((m-1)-xmean)*(transit_times[k][m]-ymean)
      denom += ((m-1)-xmean)^2
    end
    slope = numer/denom
    ttv_per[k] = slope
    intercept = ymean-slope*xmean
    for n in 1:num_times
      this_ttv = transit_times[k][n] - slope*(n-1) - intercept	#calculate ttv time relative to expected transit time
      push!(sysTTVs[k],this_ttv)
    end
    popfirst!(sysTTVs[k]) 					#removes 0.0 value from when we initialized our array
  end
  TTVFast.free_ttvfast_outputs(ttvfast_output)			#original ttv code recommends running this
  plot_ttv(transit_times,sysTTVs,ttv_per)			#plots relative ttvs over total observation time
  return sysTTVs
end

function calc_ttvfaster(sys::Any) #ExoplanetsSysSim.PlanetarySystem{ExoplanetsSysSim.Star}
#takes information from the second planet in a pair of two, and returns an array of arrays of TTVs for each planet.
  nplanets = num_planets(sys)
  sysTTVs = Array{Array{Float64,1}}(undef,nplanets)   #an array of all the TTVs for each planet
  ttv_per = zeros(nplanets)
  ttv_t0 = zeros(nplanets)
  for p in 2:nplanets
    total_transits = 1550
    jmax = 5
    mu1 = sys.planet[p-1].mass/sys.star.mass					#use reduced mass of planet and star
    mu2 = sys.planet[p].mass/sys.star.mass
    per1 = sys.orbit[p-1].P
    per2 = sys.orbit[p].P
    t01 = rand(0:per1)								#we'll use this until we have a better function for initial transit time
    t02 = rand(0:per2)
    ecc1 = sys.orbit[p-1].ecc
    ecc2 = sys.orbit[p].ecc
    omega1 = sys.orbit[p-1].omega
    omega2 = sys.orbit[p].omega
    p1 = TTVFaster.Planet_plane_hk(mu1,per1,t01,ecc1*cos(omega1),ecc1*sin(omega1)) #organize data for each planet
    p2 = TTVFaster.Planet_plane_hk(mu2,per2,t02,ecc2*cos(omega2),ecc2*sin(omega2))
    n1 = convert(Int64,floor(total_transits/per1))				#number of transits for each planet
    n2 = convert(Int64,floor(total_transits/per2))
    time1 = collect(p1.trans0 .+ range(0,stop=n1-1,length=n1) .* p1.period)
    time2 = collect(p2.trans0 .+ range(0,stop=n2-1,length=n2) .* p2.period)
    alpha0=(p1.period/p2.period)^(2//3)
    # Initialize the computation of the Laplace coefficients:
    b0=TTVFaster.LaplaceCoefficients.initialize(jmax+1,alpha0)
    # Define arrays to hold the TTVs:
    data_array = [mu1,per1,t01,ecc1*cos(omega1),ecc1*sin(omega1),mu2,per2,t02,ecc2*cos(omega2),ecc2*sin(omega2)]
    ttv_el_type = eltype(data_array) == Float64 ? Float64 : Number
    ttv1=Array{ttv_el_type}(undef,n1)
    ttv2=Array{ttv_el_type}(undef,n2)
    # Define arrays to hold the TTV coefficients and Laplace coefficients:
    f1=Array{Float64}(undef,jmax+2,5)
    f2=Array{Float64}(undef,jmax+2,5)
    b1 = Array{Float64}(undef,jmax+2,3)               #b1 was renamed from b because b is impact parameter
    TTVFaster.compute_ttv!(jmax,p1,p2,time1,time2,ttv1,ttv2,f1,f2,b1,alpha0,b0)
    #store TTVs in sysTTVs
    if p>2
      for h in 1:n1
        sysTTVs[p-1][h] += ttv1[h] #combine ttvs of planet 1 on planet 2 and planet 3 on planet 2
      end
      sysTTVs[p] = ttv2
      ttv_per[p] = per2
      ttv_t0[p] = t02
    else
      sysTTVs[1] = ttv1
      sysTTVs[2] = ttv2
      ttv_per[1] = per1
      ttv_per[2] = per2
      ttv_t0[1] = t01
      ttv_t0[2] = t02
    end
  end
  plot_ttvfaster(sysTTVs,ttv_per,ttv_t0)                         #plots data returned by calc_ttvfaster. Does not work with calc_ttvfast.
  return sysTTVs
end

function plot_ttv(transit_time::Array{Array{Float64,N} where N,1}, TTVs::Array{Array{Float64,1},1}, ttv_per::Array{Float64,1} )
#testing function that plots relative ttvs over total observation time.
  for i in 1:length(TTVs)
    this_per = round(Int64,ttv_per[i])
    plot(transit_time[i],TTVs[i], "-o", ms = 5, label = "period: $this_per")
    ylabel("TTV (days)", fontsize = 14)
    xlabel("time (days)", fontsize = 14)
    legend()
  end
end

function plot_ttvfaster(TTVs::Array{Array{Float64,1},1}, ttv_per::Array{Float64,1}, ttv_t0::Array{Float64,1})
  transit_time = Array{Array{Float64}}(undef,length(TTVs))
  for k in 1:length(TTVs)
    first_time = ttv_t0[1] + TTVs[k][1]	#we compute the first time separately so that we can overwrite the "undefined" initializer and then push the remaining times
    transit_time[k] = [first_time]
    for m in 2:length(TTVs[k][:])	#start at 2 since we already filled in the first values in the line above.
      this_time = ttv_t0[k] + ttv_per[k] * (m-1) + TTVs[k][m]	 #TODO: find out if there is a way to extract a observation-based period from ttvfaster code rather than using the period generated by SysSim
      push!(transit_time[k],this_time)
    end
  end
  plot_ttv(transit_time, TTVs, ttv_per)
end

function get_ttv_sum_stats(TTVs::Array{Array{Float64,1},1})
  num_planets = length(TTVs)
  ttv_rms = zeros(num_planets)
  ttv_iqr = zeros(num_planets)
  ttv_scatter = zeros(num_planets)
  for i in 1:num_planets
    num_events = length(TTVs[i])
    ttv_mean = mean(TTVs[i])
    abs_ttvs = zeros(num_events)
    for j in 1:num_events
      abs_ttvs[j] = abs(ttv_mean - TTVs[i][j])
    end
    ttv_scatter[i] = mean(abs_ttvs)*scatter_const
    #ttv_iqr[i] = iqr(absttvs)
    ttv_rms[i] = convert(Float64,sqrt(ttv_mean^2))
  end
  ttv_sum_stat = Dict("scatter" => ttv_scatter, "rms" => ttv_rms, "iqr" => ttv_iqr)
  #println(ttv_sum_stat)
  return ttv_sum_stat
end


text:
using Test
using AbstractPlotting
import Calculus
using GLMakie; GLMakie.activate!()
using MakieLayout
import NonlinearEquations
import Random
import SparseArrays

function callback(psi, residuals, J, i)
	if mod(i, 10) == 0
		fig, axs = PyPlot.subplots(1, 2, figsize=(16, 9))
		psi = reshape(psi, ns[2], ns[1])
		residuals = reshape(residuals, ns[2], ns[1])
		img = axs[1].imshow(map(x->x > 0 ? x : 0, psi), origin="lower", interpolation="nearest", extent=[mins[1], maxs[1], mins[2], maxs[2]])
		fig.colorbar(img, ax=axs[1])
		goodindices = filter(i->coords[2, i] == maxs[2], 1:size(coords, 2))
		img = axs[2].imshow(residuals, origin="lower", interpolation="nearest", extent=[mins[1], maxs[1], mins[2], maxs[2]])
		fig.colorbar(img, ax=axs[2])
		display(fig)
		println()
		println()
		PyPlot.close(fig)
	end
end

function plottransient(soln)
	hend = soln(tspan[end])
	scene, layout = layoutscene(resolution=(1200, 600))
	t_slider = layout[3, 1:3] = LSlider(scene; range=range(0, tspan[end]; length=500), startvalue=0)
	t_observable = lift(t->[t, t], t_slider.attributes[:value])
	h_observable = lift(t->reshape(soln(t), ns[2], ns[1])', t_slider.attributes[:value])
	residuals_observable = lift(t->reshape(f(soln(t), p, t), ns[2], ns[1])', t_slider.attributes[:value])
	ax1 = layout[1, 1] = LAxis(scene; title="Head")
	ax2 = layout[1, 2] = LAxis(scene; title="Residuals")
	xs = range(mins[1]; stop=maxs[1], length=ns[1])
	ys = range(mins[2]; stop=maxs[2], length=ns[2])
	heatmap!(ax1, xs, ys, h_observable; colorrange=extrema(hend))
	heatmap!(ax2, xs, ys, residuals_observable, colorrange=extrema(hend))
	for ax in [ax1, ax2]
		xlims!(ax, mins[1], maxs[1])
		ylims!(ax, mins[2], maxs[2])
	end
	cbar = layout[1, 3] = LColorbar(scene; limits=extrema(hend), width=30)
	timeplot = layout[2, :] = LAxis(scene; title="time", xgridvisible=false, ygridvisible=false, yticklabelsvisible=false, yticksvisible=false)
	timeplot.height=20
	lines!(timeplot, t_observable, [0, 1], linewidth=5)
	ylims!(timeplot, 0, 1)
	xlims!(timeplot, 0, tspan[end])
	display(scene)
end

function i2i1i2(i)
	i2 = 1 + mod(i - 1, ns[2])
	i1 = 1 + div(i - i2, ns[2])
	return i1, i2
end

function updateneighborarrays!(a1, a2, pair, value)
	n1, n2 = pair
	if i2i1i2(n1)[1] == i2i1i2(n2)[1]#the neighbors have the same x coordinate, update a2
		i1, i2_1 = i2i1i2(n1)
		i2_2 = i2i1i2(n2)[2]
		@assert i2_1 < i2_2
		a2[i2_1, i1] = value
	else#the neighbors have the same y coordinate, update a1
		i1_1, i2 = i2i1i2(n1)
		i1_2 = i2i1i2(n2)[1]
		@assert i1_1 < i1_2
		a1[i2, i1_1] = value
	end
end

function breakupgradient(grad)
	grad1 = zeros(ns[2], ns[1] - 1)
	grad2 = zeros(ns[2] - 1, ns[1])
	for (v, pair) in zip(grad, neighbors)
		updateneighborarrays!(grad1, grad2, pair, v)
	end
	return grad1, grad2
end

function plotgradient(grad, obsnode; vmin=minimum(grad), vmax=maximum(grad))
	grad1, grad2 = breakupgradient(grad)
	fig, axs = PyPlot.subplots(1, 2, figsize=(16, 9))
	img = axs[1].imshow(grad1, extent=[mins[1], maxs[1], mins[2], maxs[2]], origin="lower", interpolation="nearest", vmin=vmin, vmax=vmax)
	fig.colorbar(img, ax=axs[1])
	img = axs[2].imshow(grad2, extent=[mins[1], maxs[1], mins[2], maxs[2]], origin="lower", interpolation="nearest", vmin=vmin, vmax=vmax)
	fig.colorbar(img, ax=axs[2])
	display(fig)
	println()
	println()
	PyPlot.close(fig)
end

ext: jl


text:
export Operator, DerivativeMetaData, OperatorMetaData
export compute

struct Operator{𝒮, 𝒯} <: AbstractExpression
    operand::𝒮
    metadata::𝒯
end

function (o::Operator)(expr::AbstractExpression)
    return Operator(expr, o.metadata)
end

function compute(o::Operator)
    return o.metadata.operation(compute(o.operand))
end

function compute(o::Operator{𝒮, 𝒯}) where
    {𝒮 <: Nothing, 𝒯}
    return compute(o.metadata)
end

struct DerivativeMetaData{𝒪, 𝒟}
    operation::𝒪
    direction::𝒟
end

function Base.show(io::IO, o::Operator{S,T}) where
    {S <: Nothing, T <: DerivativeMetaData}
    name = Char(0x02202) * o.metadata.direction
    printstyled(io, name, color = 14 )
end

function Base.show(io::IO, o::Operator{S,T}) where
    {S <: AbstractExpression, T <: DerivativeMetaData}
    name = Char(0x02202) * o.metadata.direction
    printstyled(io, name, "(",  color = 14 )
    print(o.operand)
    printstyled(io, ")",  color = 14 )
end

struct OperatorMetaData{𝒪, 𝒩}
    operation::𝒪
    name::𝒩
end

function Base.show(io::IO, o::Operator{S,T}) where
    {S <: Nothing, T <: OperatorMetaData}
    name = o.metadata.name
    printstyled(io, name, color = 14 )
end

function Base.show(io::IO, o::Operator{S,T}) where
    {S <: AbstractExpression, T <: OperatorMetaData}
    name = o.metadata.name
    printstyled(io, name, "(",  color = 14 )
    print(o.operand)
    printstyled(io, ")",  color = 14 )
end

to_expr(t::Operator) = Expr(:call, t, to_expr(t.operand))

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
"""
Author: David Rushing Dewhurst

Copyright: David Rushing Dewhurst and Charles River Analytics Inc., 2021 - present

Released under the MIT license.
"""

module CrimsonSkyline

using DataFrames
using DataStructures
using Distributions
using PrettyPrint: pprintln
using StatsBase
using Logging
using JuliaDB
using JSON
using MacroTools
using MacroTools: prewalk
using Plots
using PDMats
using SQLite

include("modeling/util.jl")

include("modeling/trace.jl")
include("modeling/distributions.jl")

include("representation/results.jl")

include("modeling/field.jl")

include("representation/statistics.jl")

include("modeling/effects.jl")

include("inference/forward.jl")
include("inference/rejection.jl")
include("inference/importance.jl")
include("inference/metropolis.jl")
include("inference/nested.jl")

include("representation/io.jl")
include("representation/plot.jl")

include("inference/kernel.jl")
include("inference/inference.jl")
include("representation/db.jl")

include("modeling/dsl.jl")

end # module


text:
using CUDA

a = CUDA.DEVICE_ATTRIBUTE_MULTI_GPU_BOARD
@show multigpuboard = CUDA.attribute(device(), a)

a = CUDA.DEVICE_ATTRIBUTE_MEMORY_CLOCK_RATE
@show memclockrate = CUDA.attribute(device(), a)

a = CUDA.DEVICE_ATTRIBUTE_MAX_SHARED_MEMORY_PER_BLOCK
@show shmemperblock = CUDA.attribute(device(), a)

a = CUDA.CU_DEVICE_ATTRIBUTE_MAX_SHARED_MEMORY_PER_MULTIPROCESSOR
@show shmempersm = CUDA.attribute(device(), a)

a = CUDA.CU_DEVICE_ATTRIBUTE_MAX_REGISTERS_PER_BLOCK
@show regsperblock = CUDA.attribute(device(), a)

a = CUDA.CU_DEVICE_ATTRIBUTE_MAX_THREADS_PER_MULTIPROCESSOR
@show threadspersm = CUDA.attribute(device(), a)

a = CUDA.CU_DEVICE_ATTRIBUTE_CLOCK_RATE
CUDA.attribute(device(), a)

CUDA.CU_DEVICE_ATTRIBUTE_MEMORY_CLOCK_RATE
CUDA.attribute(device(), a)

CUDA.CU_DEVICE_ATTRIBUTE_MAX_THREADS_PER_MULTIPROCESSOR
CUDA.attribute(device(), a)

CUDA.CU_DEVICE_ATTRIBUTE

a = 0
CUDA.cuDeviceGetProperties(a, device())

b = CUDA.CUdevprop_st


text:
#*****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

using TensorFlow
import TensorFlow.API: shape, gather, concat, expand_dims, div_

type Normalizer{T<:AbstractFloat}
    ndim::Int64
    maxes::Array{T}
    mins::Array{T}
end

function Normalizer(data::TFDataset; converttype::Union{Type,Void}=nothing)
    N = ndims(data.X)
    rmdims = collect(1:N-1)
    #find max/mins per column over entire dataset
    maxes = maximum(data.X, rmdims)
    mins = minimum(data.X, rmdims)
    if converttype != nothing
        maxes = convert(Array{converttype}, maxes)
        mins = convert(Array{converttype}, mins)
    end
    norm = Normalizer(N, maxes, mins)
    norm
end

#"""
#"""
function normalize01(norm::Normalizer, Xin::Tensor)
    N = norm.ndim
    maxes = constant(norm.maxes)
    mins = constant(norm.mins)
    tileshape = concat(Tensor(0), Tensor([gather(shape(Xin), Tensor(collect(0:N-2))); Tensor([1])]))
    tiledmaxes = tile(maxes, tileshape)
    tiledmins = tile(mins, tileshape)
    Xout = div_((Xin - tiledmins), (tiledmaxes - tiledmins))
    Xout
end

function normalize(norm::Normalizer, Xin::Tensor, minval::AbstractFloat=-1.0,
    maxval::AbstractFloat=1.0)

    X01 = normalize(Xin)
    max = constant(maxval)
    min = constant(minval)
    shape = get_shape(Xin)
    tiledmax = tile(max, Tensor(shape))
    tiledmin = tile(min, Tensor(shape))
    Xout = X01 * (tiledmax - tiledmin) + tiledmin
    Xout
end


text:
# cudnn wrapper

"""

`conv_f(x,w;kwargs...)` executes convolutions or cross-correlations
using filters specified with `w` over tensor `x`.
Here is a description of all available keyword arguments:
* padding: the number of extra zeros implicitly concatenated at the start and at the end of each dimension. Default=floor((filterSize-1)/2) which preserves the input size when filterSize is odd and stride=1.
* stride: the number of elements to slide to reach the next filtering window. Default=1.
* upscale: upscale factor for each dimension. Default=1.
* mode: 0 for convolution and 1 for cross-correlation.  Default=0.
* alpha: can be used to scale the result. Default=1.
* algo: specifies which convolution algorithm shoud be used to compute the results. Default=0. See the CUDNN User Guide for details.
* workSpace: data pointer to GPU memory to a workspace needed to able to execute the specified algorithm. Default=C_NULL.
* workSpaceSizeInBytes: the size in bytes of the provided workSpace. Default=0.
* handle: handle to a previously created cuDNN context. Default=Dex allocated context.

"""
function conv_f{T}(x::DexArray{T},w::DexArray{T};
                  handle=cudnnhandle, alpha=one(T), beta=zero(T),
                  algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0, o...)
    y = similar(x, cdims(x,w;o...))
    @cuda(cudnn, cudnnConvolutionForward,
          (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,UInt32,Cptr,Csize_t,Ptr{T},Cptr,Ptr{T}),
          handle,Ref(alpha),TD(x),x,FD(w),w,CD(x,w;o...),algo,workSpace,workSpaceSizeInBytes,Ref(beta),TD(y),y)
    return y
end

"""

`conv_bx(x,w,dy;kwargs...)` This function computes the convolution gradient with
respect to the output tensor, returning results in gradDesc.

"""
function conv_bx{T}(x::DexArray{T},w::DexArray{T},dy::DexArray{T};
                   handle=cudnnhandle, alpha=one(T), beta=zero(T),
                   algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0, o...)
    dx = similar(x)
    @cuda(cudnn,cudnnConvolutionBackwardData,
          (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
          handle,Ref(alpha),FD(w),w,TD(dy),dy,CD(x,w;o...),algo,workSpace,workSpaceSizeInBytes,Ref(beta),TD(dx),dx)
    return dx
end

"""

`conv_bw(x,w,dy;kwargs...)` This function computes the convolution gradient with
 respect to filter coefficients, returning results in gradDesc.

"""
function conv_bw{T}(x::DexArray{T},w::DexArray{T},dy::DexArray{T};
                   handle=cudnnhandle, alpha=one(T), beta=zero(T),
                   algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0, o...)
    dw = similar(w)
    @cuda(cudnn,cudnnConvolutionBackwardFilter,
          (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,     UInt32,Cptr,     Csize_t,             Ptr{T},Cptr,Ptr{T}),
          handle,Ref(alpha),TD(x),x,TD(dy),dy,CD(x,w;o...),algo,workSpace,workSpaceSizeInBytes,Ref(beta),FD(dw),dw)
    return dw
end

"""

`pool_f(x;kwargs...)` computes pooling of input values (i.e., the
maximum or average of several adjacent values) to produce an output
with smaller height and/or width.

Here is a description of all available keyword arguments:
* window: the pooling window size for each dimension. Default=2.
* padding: the number of extra zeros implicitly concatenated at the start and at the end of each dimension. Default=0.
* stride: the number of elements to slide to reach the next pooling window. Default=same as window.
* mode: 0 for max, 1 for average including padded values, 2 for average excluding padded values.  Default=0.
* maxpoolingNanOpt: Nan numbers are not propagated if 0, they are propagated if 1. Default=0.
* alpha: can be used to scale the result. Default=1.
* handle: Handle to a previously created cuDNN context. Default=Dex allocated context.

"""

function pool_f{T}(x::DexArray{T}; handle=cudnnhandle, alpha=one(T), beta=zero(T), o...)
    y = similar(x, pdims(x; o...))
    @cuda(cudnn, cudnnPoolingForward,
          (Cptr, Cptr,      Ptr{T},    Cptr,Ptr{T},Ptr{T},   Cptr,Ptr{T}),
          handle,PD(x;o...),Ref(alpha),TD(x),x,    Ref(beta),TD(y),y)
    return y
end

function pool_b{T}(x::DexArray{T},y::DexArray{T},dy::DexArray{T};
                  handle=cudnnhandle, alpha=one(T), beta=zero(T), o...)
    dx = similar(x)
    @cuda(cudnn,cudnnPoolingBackward,
          (Cptr,Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Ptr{T},Cptr,Ptr{T}),
          handle,PD(x;o...),Ref(alpha),TD(y),y,TD(dy),dy,TD(x),x,Ref(beta),TD(dx),dx)
    return dx
end

"""

`dropout_f(x,dropout;kwargs...)` computes the forward dropout operation over input values
to produce an output.

Here is a description of all available keyword arguments:
* handle: Handle to a previously created cuDNN context. Default=Dex allocated context.
* dropout: The probability with which the value from input would be propagated through the dropout layer.
* reserveSpace: Pointer to user-allocated GPU memory used by this function. It is expected that contents of reserveSpace doe not change between cudnnDropoutForward and cudnnDropoutBackward calls.
* reserveSpaceSizeInBytes: Specifies size in bytes of the provided memory for the reserve space
TODO: find a neat way to hide reserveSpace. through some kind of global variable managerment.

"""

function dropout_f{T}(x::DexArray{T};dropout=0.5,handle=cudnnhandle)
    y = similar(x)
    reservesize_p = Csize_t[0]
    @cuda(cudnn, cudnnDropoutGetReserveSpaceSize,(Cptr,Ptr{Csize_t}),TD(x),reservesize_p)
    reservesize = reservesize_p[1]
    reservespace = DexArray(Int8, Int(reservesize))
    @cuda(cudnn, cudnnDropoutForward,
          (Cptr,Cptr,        Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Csize_t),
          handle,DD(dropout=dropout),TD(x),x,   TD(y),y,    reservespace,reservesize)
    return y, reservespace  # reservespace should be hide...
end

function dropout_b{T}(dy::DexArray{T},reservespace,dropout::Float64;handle=cudnnhandle)
    dx = similar(dy)
    @cuda(cudnn, cudnnDropoutBackward,
          (Cptr,Cptr,                Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Csize_t),
          handle,DD(dropout=dropout),TD(dy),dy, TD(dx),dx,  reservespace,length(reservespace))
    return dx
end

"""

`rnn_ft(x,w,seqlength;kwargs...)` executes the forward training of the recurrent neural network and produce an output y.
* seqlength: Number of iterations to unroll over.
Here is a description of all available keyword arguments:
* handle: Handle to a previously created cuDNN context. Default=Dex allocated context.
* hiddensize: size of the hidden unit

"""

function rnn_ft{T}(x::DexArray{T};handle=cudnnhandle,hiddensize=128,numlayers=1,mode="relu",bidirectional=0,dropout=0.5)
    seqlength = size(x,1)
    batchsize = size(x,2)
    hidden_shape = (numlayers,size(x,2),(bidirectional==1?2:1)*hiddensize)
    info("hidden_shape: $hidden_shape")
    hx = C_NULL # DexArray(rand(T,hidden_shape)) # use rand to initial hidden state of the RNN
    cx = C_NULL # DexArray(rand(T,hidden_shape)) # use rand to initial cell state of the network
    hy = C_NULL # similar(hx)
    cy = C_NULL # similar(cx)
    y = DexArray(T,seqlength,batchsize,(bidirectional == 1?2:1)*hiddensize)
    ydescs = fill(TD(DexArray(T,size(y,2),size(y,3),1)),seqlength)
    xdesc = TD(DexArray(T,size(x,2),size(x,3),1))
    xdescs = fill(xdesc,seqlength)
    rnndesc = RD(x,hiddensize,numlayers,mode,bidirectional,dropout)
    worksize_p = Cptr[0]
    @cuda(cudnn, cudnnGetRNNWorkspaceSize,(Cptr,Cptr,Cint,Ptr{Cptr},Cptr),handle,rnndesc,Cint(seqlength),xdescs,worksize_p)
    worksize = worksize_p[1]
    info("worksize: $(Int(worksize))")
    workspace = DexArray(Int8, Int(worksize))
    resevesize_p = Cptr[0]
    @cuda(cudnn, cudnnGetRNNTrainingReserveSize,(Cptr,Cptr,Cint,Ptr{Cptr},Cptr),handle,rnndesc,Cint(seqlength),xdescs,resevesize_p)
    reservesize = resevesize_p[1]
    info("reservesizesize: $(Int(reservesize))")
    reservespace = DexArray(Int8, Int(reservesize))
    wsize_p = Cint[0]
    @cuda(cudnn, cudnnGetRNNParamsSize,(Cptr,Cptr,Cptr,Ptr{Cint},UInt32),handle,rnndesc,xdesc,wsize_p,DT(x))
    wsize = wsize_p[1]
    info("wsize: $(Int(wsize/T.size))")
    w = DexArray(rand(T,Int(wsize/T.size))) # w should be 1d array
    wdesc = FD(DexArray(T,Int(wsize/T.size),1,1))
    @cuda(cudnn, cudnnRNNForwardTraining,
          (Cptr,Cptr,Cint,Ptr{Cptr},Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},
          Ptr{Cptr},Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Cint,Cptr,Cint),
          handle,rnndesc,Cint(seqlength),xdescs,x,C_NULL,hx,C_NULL,cx,wdesc,w,
          ydescs,y,C_NULL,hy,C_NULL,cy,workspace,Cint(worksize),reservespace,Cint(reservesize))
    # @cuda(cudnn, cudnnRNNForwardTraining,
    #       (Cptr,Cptr,Cint,Ptr{Cptr},Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},
    #       Ptr{Cptr},Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Cint,Cptr,Cint),
    #       handle,rnndesc,Cint(seqlength),xdescs,x,TD(hx),hx,TD(cx),cx,wdesc,w,
    #       ydescs,y,TD(hy),hy,TD(cy),cy,workspace,Cint(worksize),reservespace,Cint(reservesize))
    return y, hy, cy, w, rnndesc, reservespace
end

"""

`rnn_fi(rnndesc,x,w;kwargs...)` executes the forward inference of the recurrent neural network and produce an output y.

"""

function rnn_fi{T}(rnndesc,x::DexArray{T},w::DexArray{T},seqlength::Int;handle=cudnnhandle,hiddensize=128,numlayers=1,mode="relu",bidirectional=0,dropout=0)
    seqlength = size(x,1)

    y = DexArray(T,seqlength,batchsize,(bidirectional == 1)?hiddensize*2:hiddensize)
    # hidden shape: number of layers, batchsize(2nd dim of x), hiddensize(*2 if bidirectional)
    hidden_shape = (numlayers,size(x,2),(bidirectional == 1)?hiddensize*2:hiddensize)
    hx = DexArray(rand(hidden_shape)) # use rand to initial hidden state of the RNN
    cx = DexArray(rand(hidden_shape)) # use rand to initial cell state of the network
    hy = similar(hx)
    cy = similar(cx)

    worksize_p = Cptr[0]
    @cuda(cudnn, cudnnGetRNNWorkspaceSize,(Ptr{Cptr},Cptr),xdescs,worksize_p)
    worksize = worksize_p[1]
    workspace = DexArray(Int8, Int(worksize))

    @cuda(cudnn, cudnnRNNForwardInference,
          (Cptr,Cptr,Cint,                                                     Ptr{Cptr},Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Csize_t),
          handle,rnndesc,Cint(seqlength),xdescs,x,        TD(hx),hx,  TD(cx),cx,  FD(w),w,    TD(y),y,    TD(hy),hy,  TD(cy),cy,  workspace,worksize)
    return y, hy, cy
end

"""
`rnn_bx(x;kwargs...)` compute the geadients dy, dhy, dhc, weights w and input gradients dx, dhx, dcx.
workspace is required for intermediate storage. The data in reserveSpace must have previously been generated by cudnnRNNForwardTraining.

"""

function rnn_bx{T}(rnndesc,x::DexArray{T},y::DexArray{T},w::DexArray{T},dy::DexArray{T},dhy::DexArray{T},dcy::DexArray{T},reservespace,seqlength::Int;handle=cudnnhandle,hiddensize=128,numlayers=1,mode="relu",bidirectional=0,dropout=0)
    ydescs = fill(TD(y),seqlength)
    dydescs = fill(TD(dy),seqlength)
    dxdescs = fill(TD(x),seqlength)
    hx = DexArray(rand(size(dhy))) # use rand to initial hidden state of the RNN
    cx = DexArray(rand(size(dcy))) # use rand to initial cell state of the network

    dx = similar(x)
    dhx = similar(hx)
    dcx = similar(cx)
    worksize_p = Cptr[0]
    @cuda(cudnn, cudnnGetRNNWorkspaceSize,(Ptr{Cptr},Cptr),xdescs,worksize_p)
    worksize = worksize_p[1]
    workspace = DexArray(Int8, Int(worksize))
    @cuda(cudnn, cudnnRNNBackwardData,
          (Cptr,Cptr,Cint,Ptr{Cptr},Ptr{T},    Ptr{Cptr},Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Ptr{Cptr},Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Csize_t,Cptr,Csize_t),
          handle,rnndesc,Cint(seqlength),ydescs,y,dydescs,dy,      TD(dhy),dhy,TD(dcy),dcy,FD(w),w,    TD(hx),hx,  TD(cx),cx,  dxdescs,dx,      TD(dhx),dhx,TD(dcx),dcx,workspace,worksize,reservespace,length(reservespace))
    return dx, dhx, dxc, reservespace
end

"""

`rnn_bw(x;kwargs...)` accumulates weight gradients dw from the recurrent neural network

"""

function rnn_bw{T}(x::DexArray{T},y::DexArray{T},w::DexArray{T},dw::DexArray{T},reservespace,seqlength::Int;handle=cudnnhandle,hiddensize=128,numlayers=1,mode="relu",bidirectional=0)
    xdescs = fill(TD(x),seqlength)
    ydescs = fill(TD(y),seqlength)

    worksize_p = Cptr[0]
    @cuda(cudnn, cudnnGetRNNWorkspaceSize,(Ptr{Cptr},Cptr),xdescs,worksize_p)
    worksize = worksize_p[1]
    workspace = DexArray(Int8, Int(worksize))
    @cuda(cudnn, cudnnRNNBackwardWeights,
          (Cptr,Cint,Ptr{Cptr},Ptr{T}    ,Cptr,Ptr{T},Ptr{Cptr},Ptr{T},Cptr,Csize_t,Cptr,Ptr{T}    ,Cptr,Csize_t),
          handle,Cint(seqlength),xdescs,x,TD(hx),hx,  ydescs,y,        workspace,worksize,FD(dw),dw,reservespace,length(reservespace))
    return dw
end

# cudnn datatype
# an enumerated type indicating the data type to which a tensor descriptor
# or filter descriptor refers.
DT(::DexArray{Float32})=UInt32(0)
DT(::DexArray{Float64})=UInt32(1)
DT(::DexArray{Float16})=UInt32(2)

# accquire dimension
function cdims{T,N}(x::DexArray{T,N},w::DexArray{T,N}; padding=0, stride=1, o...)
    ntuple(N) do i
        if i < N-1
            pi = if isa(padding,Number); padding; else padding[i]; end
            si = if isa(stride,Number); stride; else stride[i]; end
            1 + div(size(x,i) - size(w,i) + 2*pi, si)
        elseif i == N-1
            size(w,N)
        else # i == N
            size(x,N)
        end
    end
end

function pdims{T,N}(x::DexArray{T,N}; window=2, padding=0, stride=window, o...)
    ntuple(N) do i
        if i < N-1
            wi = (if isa(window,Number); window; else window[i]; end)
            pi = (if isa(padding,Number); padding; else padding[i]; end)
            si = (if isa(stride,Number); stride; else stride[i]; end)
            1 + div(size(x,i) + 2*pi - wi, si)
        else
            size(x,i)
        end
    end
end

# cudnn descriptors

type TD; ptr
    function TD(a::DexArray)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreateTensorDescriptor,(Ptr{Cptr},),d)
        n = ndims(a)
        sz = [Cint(size(a,n-i+1)) for i=1:n]    # why dim reversed?
        st = [Cint(stride(a,n-i+1)) for i=1:n]
        @cuda(cudnn,cudnnSetTensorNdDescriptor,
              (Cptr,UInt32,Cint,Ptr{Cint},Ptr{Cint}),
              d[1], DT(a), n, sz, st)
        td = new(d[1])
        finalizer(td, x->@cuda(cudnn,cudnnDestroyTensorDescriptor,(Cptr,),x.ptr))
        return td
    end
end

type FD; ptr
    function FD(a::DexArray)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreateFilterDescriptor,(Ptr{Cptr},),d)
        n = ndims(a)
        sz = [Cint(size(a,n-i+1)) for i=1:n]
        @cuda(cudnn,cudnnSetFilterNdDescriptor,
              (Cptr,UInt32,UInt32,Cint,Ptr{Cint}),
              d[1], DT(a), 0,     n,   sz)
        fd = new(d[1])
        finalizer(fd, x->@cuda(cudnn,cudnnDestroyFilterDescriptor,(Cptr,),x.ptr))
        return fd
    end
end

type CD; ptr
    function CD(x::DexArray,w::DexArray; padding=0, stride=1, upscale=1, mode=0)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreateConvolutionDescriptor,(Ptr{Cptr},),d)
        nd = ndims(x)-2
        @cuda(cudnn,cudnnSetConvolutionNdDescriptor,
              (Cptr,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},UInt32,UInt32),
              d[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(upscale,nd),mode,DT(x))
        cd = new(d[1])
        finalizer(cd, x->@cuda(cudnn,cudnnDestroyConvolutionDescriptor,(Cptr,),x.ptr))
        return cd
    end
end

type PD; ptr
    function PD(x::DexArray; window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=0)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreatePoolingDescriptor,(Ptr{Cptr},),d)
        nd = ndims(x)-2
        @cuda(cudnn,cudnnSetPoolingNdDescriptor,
              (Cptr,UInt32,UInt32,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint}),
              d[1],mode,maxpoolingNanOpt,nd,cdsize(window,nd),cdsize(padding,nd),cdsize(stride,nd))
        pd = new(d[1])
        finalizer(pd, x->@cuda(cudnn,cudnnDestroyPoolingDescriptor,(Cptr,),x.ptr))
        return pd
    end
end

type DD; ptr
    function DD(;dropout=0.5,handle=cudnnhandle)
        d = Cptr[0]
        @cuda(cudnn, cudnnCreateDropoutDescriptor, (Ptr{Cptr},),d)
        statessize_p = Csize_t[0]
        @cuda(cudnn, cudnnDropoutGetStatesSize, (Cptr,Ptr{Csize_t}), handle,statessize_p)
        statessize = statessize_p[1]
        states = DexArray(Int8, Int(statessize))
        @cuda(cudnn, cudnnSetDropoutDescriptor,
              (Cptr,Cptr, Cfloat,Cptr,  Cint,Culonglong),
              d[1],handle,Cfloat(dropout),states,statessize, 0)
        dd = new(d[1])
        finalizer(dd, x->@cuda(cudnn, cudnnDestroyDropoutDescriptor, (Cptr,),x.ptr))
        return dd
    end
end

rnnmode = Dict("relu"=>0,"tanh"=>1,"lstm"=>2,"GRU"=>3)
type RD; ptr
    function RD(x::DexArray,hiddensize::Int,numlayers::Int,mode::AbstractString,bidirectional::Int,dropout::Float64)
        d = Cptr[0]
        @cuda(cudnn, cudnnCreateRNNDescriptor,(Ptr{Cptr},),d)

        @cuda(cudnn, cudnnSetRNNDescriptor,
              (Cptr,Cint, Cint,                     Cptr,Cint,Cint,Cint,UInt32),
              d[1],Cint(hiddensize),Cint(numlayers),DD(dropout=dropout),0,Cint(bidirectional),Cint(rnnmode[mode]),DT(x))
        rd = new(d[1])
        finalizer(rd, x->@cuda(cudnn, cudnnDestroyRNNDescriptor, (Cptr,),x.ptr))
        return rd
    end
end

import Base: unsafe_convert
unsafe_convert(::Type{Cptr}, td::TD)=td.ptr
unsafe_convert(::Type{Cptr}, fd::FD)=fd.ptr
unsafe_convert(::Type{Cptr}, cd::CD)=cd.ptr
unsafe_convert(::Type{Cptr}, pd::PD)=pd.ptr
unsafe_convert(::Type{Cptr}, dd::DD)=dd.ptr
unsafe_convert(::Type{Cptr}, rd::RD)=rd.ptr

function cdsize(w, nd)
    if isa(w,Integer)
        fill(Cint(w),nd)
    elseif length(w)==nd
        [ Cint(w[nd-i+1]) for i=1:nd ]
    else
        throw(DimensionMismatch())
    end
end

# convolution padding size that preserves the input size when filter size is odd and stride=1
padsize(w)=ntuple(i->div(size(w,i)-1,2), ndims(w)-2)


text:
using Base: UUID
using Preferences, Test, TOML

function activate(f::Function, env_dir::String)
    saved_active_project = Base.ACTIVE_PROJECT[]
    Base.ACTIVE_PROJECT[] = env_dir
    try
        f()
    finally
        Base.ACTIVE_PROJECT[] = saved_active_project
    end
end

function with_temp_depot(f::Function)
    mktempdir() do dir
        saved_depot_path = copy(Base.DEPOT_PATH)
        empty!(Base.DEPOT_PATH)
        push!(Base.DEPOT_PATH, dir)
        try
            f()
        finally
            empty!(Base.DEPOT_PATH)
            append!(Base.DEPOT_PATH, saved_depot_path)
        end
    end
end

function activate_and_run(project_dir::String, code::String; env::Dict = Dict())
    mktempdir() do dir
        open(joinpath(dir, "test_code.jl"), "w") do io
            write(io, code)
        end

        out = Pipe()
        cmd = setenv(`$(Base.julia_cmd()) --project=$(project_dir) $(dir)/test_code.jl`,
                     env..., "JULIA_DEPOT_PATH" => Base.DEPOT_PATH[1])
        p = run(pipeline(cmd, stdout=out, stderr=out); wait=false)
        close(out.in)
        wait(p)
        output = String(read(out))
        if !success(p)
            println(output)
        end
        @test success(p)
        return output
    end
end

# Some useful constants
up_uuid = UUID(TOML.parsefile(joinpath(@__DIR__, "UsesPreferences", "Project.toml"))["uuid"])
up_path = joinpath(@__DIR__, "UsesPreferences")

@testset "Preferences" begin
    # Ensure there is no LocalPreferences.toml file in UsesPreferences:
    local_prefs_toml = joinpath(up_path, "LocalPreferences.toml")
    rm(local_prefs_toml; force=true)
    with_temp_depot() do
        # Start with the default test of the backend being un-set, we just get the default
        activate_and_run(up_path, """
            using UsesPreferences, Test, Preferences
            using Base: UUID
            @test load_preference($(repr(up_uuid)), "backend") === nothing
            @test UsesPreferences.backend == "OpenCL"
        """)

        # Next, change a setting
        activate_and_run(up_path, """
            using UsesPreferences
            UsesPreferences.set_backend("CUDA")
        """)

        # Ensure that's showing up in LocalPreferences.toml:
        prefs = TOML.parsefile(local_prefs_toml)
        @test haskey(prefs, "UsesPreferences")
        @test prefs["UsesPreferences"]["backend"] == "CUDA"

        # Now show that it forces recompilation
        did_precompile(output) = occursin("Precompiling UsesPreferences [$(string(up_uuid))]", output)
        cuda_test = """
        using UsesPreferences, Test
        @test UsesPreferences.backend == "CUDA"
        """
        output = activate_and_run(up_path, cuda_test; env=Dict("JULIA_DEBUG" => "loading"))
        @test did_precompile(output)

        # Show that it does not force a recompile the second time
        output = activate_and_run(up_path, cuda_test; env=Dict("JULIA_DEBUG" => "loading"))
        @test !did_precompile(output)

        # Test non-compiletime preferences a bit
        activate_and_run(up_path, """
            using UsesPreferences, Test, Preferences
            using Base: UUID
            @test load_preference($(repr(up_uuid)), "username") === nothing
            @test !UsesPreferences.has_username()
            @test UsesPreferences.get_username() === nothing
            UsesPreferences.set_username("giordano")
            @test UsesPreferences.get_username() == "giordano"
        """)

        # This does not cause a recompilation, and we can also get the username back again:
        username_test = """
        using UsesPreferences, Test, Preferences
        @test UsesPreferences.get_username() == "giordano"
        """
        output = activate_and_run(up_path, username_test; env=Dict("JULIA_DEBUG" => "loading"))
        @test !did_precompile(output)

        _prefs_delete_username = "delete_preferences!($(repr(up_uuid)), \"username\"; block_inheritance = false)"
        _delete_username = "UsesPreferences.delete_username()"
        _has_username = "@test UsesPreferences.has_username()"
        _doesnt_have_username = "!@test UsesPreferences.has_username()"
        _set_username = "UsesPreferences.set_username(\"giordano\")"
        _get_username = "@test UsesPreferences.get_username() == \"giordano\""
        _doesnt_have_set_get_username = """
        $(_doesnt_have_username)
        $(_set_username)
        $(_get_username)
        """
        snippets = [
            _prefs_delete_username,
            _doesnt_have_set_get_username,
            _has_username,
            _delete_username,
            _doesnt_have_set_get_username,
            _has_username,
            _prefs_delete_username,
            _doesnt_have_set_get_username,
            _has_username,
            _prefs_delete_username,
            _doesnt_have_set_get_username,
            _has_username,
        ]
        for snippet in snippets
            code = """
            using UsesPreferences, Test, Preferences
            using Base: UUID

            $(snippet)
            """
            output = activate_and_run(up_path, username_test; env=Dict("JULIA_DEBUG" => "loading"))
            @test !did_precompile(output)
        end
    end
end

# Load UsesPreferences, as we need it loaded for some set/get trickery below
activate(up_path) do
    eval(:(using UsesPreferences))
end
@testset "Inheritance" begin
    # Ensure there is no LocalPreferences.toml file in UsesPreferences:
    local_prefs_toml = joinpath(up_path, "LocalPreferences.toml")
    rm(local_prefs_toml; force=true)
    with_temp_depot() do
        mktempdir() do env_dir
            # We're going to create a higher environment
            push!(Base.LOAD_PATH, env_dir)

            # We're going to add `UsesPreferences` to this environment
            open(joinpath(env_dir, "Project.toml"), "w") do io
                TOML.print(io, Dict(
                    "deps" => Dict(
                        "UsesPreferences" => string(up_uuid),
                    )
                ))
            end

            # We're going to write out some Preferences for UP in the higher environment
            activate(env_dir) do
                set_preferences!(up_uuid, "location" => "outer_public"; export_prefs=true)
                # Verify that this is stored in the environment's Project.toml file
                proj = Base.parsed_toml(joinpath(env_dir, "Project.toml"))
                @test haskey(proj, "preferences")
                @test haskey(proj["preferences"], "UsesPreferences")
                @test proj["preferences"]["UsesPreferences"]["location"] == "outer_public"
                @test load_preference(up_uuid, "location") == "outer_public"

                # Add preferences to the outer env's `LocalPreferences.toml`
                set_preferences!(up_uuid, "location" => "outer_local")
                prefs = Base.parsed_toml(joinpath(env_dir, "LocalPreferences.toml"))
                @test haskey(prefs, "UsesPreferences")
                @test prefs["UsesPreferences"]["location"] == "outer_local"
                @test load_preference(up_uuid, "location") == "outer_local"
            end

            # Ensure that we can load the preferences the same even if we exit the `activate()`
            @test load_preference(up_uuid, "location") == "outer_local"

            # Next, we're going to create a lower environment, add some preferences there, and ensure
            # the inheritance works properly.
            activate(up_path) do
                # Ensure that activating this other path doesn't change anything
                @test load_preference(up_uuid, "location") == "outer_local"

                # Set a local preference in this location, which should be the first location on the load path
                set_preferences!(up_uuid, "location" => "inner_local")
                prefs = Base.parsed_toml(joinpath(up_path, "LocalPreferences.toml"))
                @test haskey(prefs, "UsesPreferences")
                @test prefs["UsesPreferences"]["location"] == "inner_local"
                @test load_preference(up_uuid, "location") == "inner_local"
            end

            # Let's add some complex objects, test that recursive merging works, and that
            # the special meaning of `nothing` and `missing` works
            activate(env_dir) do
                set_preferences!(up_uuid, "nested" => Dict(
                    "nested2" => Dict("a" => 1, "b" => 2),
                    "leaf" => "hello",
                ); export_prefs=true)
                set_preferences!(up_uuid, "nested" => Dict(
                    "nested2" => Dict("b" => 3)),
                )

                nested = load_preference(up_uuid, "nested")
                @test isa(nested, Dict) && haskey(nested, "nested2")
                @test nested["nested2"]["a"] == 1
                @test nested["nested2"]["b"] == 3
                @test nested["leaf"] == "hello"
            end

            # Add another layer in the inner environment
            activate(up_path) do
                set_preferences!(up_uuid, "nested" => Dict(
                    "nested2" => Dict("a" => "foo"),
                    "leaf" => "world",
                ))
                nested = load_preference(up_uuid, "nested")
                @test isa(nested, Dict) && haskey(nested, "nested2")
                @test nested["nested2"]["a"] == "foo"
                @test nested["nested2"]["b"] == 3
                @test nested["leaf"] == "world"
            end

            # Set the local setting of the upper environment to `missing`; this causes it to
            # pass through and `b` will suddenly equal `2`:
            activate(env_dir) do
                # Test that trying to over-set a preference in another package fails unless we force it
                @test_throws ArgumentError set_preferences!(up_uuid, "nested" => nothing)
                set_preferences!(up_uuid, "nested" => Dict(
                    "nested2" => missing,
                    "leaf" => nothing,
                ); force=true)
                nested = load_preference(up_uuid, "nested")
                @test isa(nested, Dict) && haskey(nested, "nested2")
                @test nested["nested2"]["a"] == 1
                @test nested["nested2"]["b"] == 2

                # Let's check that the `__clear__` keys are what we expect:
                prefs = Base.parsed_toml(joinpath(env_dir, "LocalPreferences.toml"))
                @test prefs["UsesPreferences"]["nested"]["__clear__"] == ["leaf"]
                @test !haskey(prefs["UsesPreferences"]["nested"], "leaf")
                @test !haskey(nested, "leaf")
            end

            # Show that it cascades down to the lower levels as well
            activate(up_path) do
                nested = load_preference(up_uuid, "nested")
                @test isa(nested, Dict) && haskey(nested, "nested2")
                @test nested["nested2"]["a"] == "foo"
                @test nested["nested2"]["b"] == 2
                @test nested["leaf"] == "world"
            end
        end
    end
end


text:
struct Texture
    """
      0 -- means texture is not embedded
    ≥ 1 -- texture is embedded and located in [scene.mTextures] array
    """
    type::aiTextureType
    embedded::Int64
    file::String
    data::Array
end

@with_kw struct Material
    ambient::aiColor4D = aiColor4D(0, 0, 0, 0)
    diffuse::aiColor4D = aiColor4D(0, 0, 0, 0)
    specular::aiColor4D = aiColor4D(0, 0, 0, 0)
    emissive::aiColor4D = aiColor4D(0, 0, 0, 0)

    textures::Dict{aiTextureType, Vector{Texture}} # texture type => [Texture]

    shininess::Float32 = 0f0
    shinpercent::Float32 = 0f0
    opacity::Float32 = 0f0
    transparencyfactor::Float32 = 0f0
    reflectivity::Float32 = 0f0
    refracti::Float32 = 0f0
    bumpscaling::Float32 = 0f0
end

struct MetaMesh
    mesh::GeometryBasics.Mesh
    material_idx::Int64
    tangents::Union{Vector{Vec3f0}, Missing}
    bitangents::Union{Vector{Vec3f0}, Missing}
end

struct Node
    meshes::Vector{MetaMesh}
    children::Vector{Node}
end

struct Scene
    node::Node
    materials::Dict{UInt32, Material}
    path::String
end


text:
function getnodestopology(nodes, faces)
    for node in nodes
        for face in faces
            if node.index in face.inodes
                push!(node.ifaces, face.index)
                if !(face.iowner in node.icells)
                    push!(node.icells, face.iowner)
                end
            end
        end
    end
end

function getfacestopology(nodes, faces)
end

function getcellstopology(nodes, faces, cells)
    for i in 1:length(faces)
        findnodes = false
        iowner = faces[i].iowner
        ineighbour = faces[i].ineighbour
        if !(i in cells[iowner].ifaces)
            push!(cells[iowner].ifaces, i)
            findnodes = true
            for node in faces[i].inodes
                if !(node in cells[iowner].inodes)
                    push!(cells[iowner].inodes, node)
                end
            end
        end
        if ineighbour != -1
            push!(cells[iowner].ineighbours, ineighbour)
            push!(cells[ineighbour].ineighbours, iowner)
            if !(i in cells[ineighbour].ifaces)
                push!(cells[ineighbour].ifaces, i)
                findnodes = true
                for node in faces[i].inodes
                    if !(node in cells[ineighbour].inodes)
                        push!(cells[ineighbour].inodes, node)
                    end
                end
            end
        end
    end
    for cell in cells
        cell.nneighbours = length(cell.ineighbours)
        for iface in cell.ifaces
            if faces[iface].iowner == cell.index
                push!(cell.facesign, 1)
            elseif faces[iface].ineighbour == cell.index
                push!(cell.facesign, -1)
            end
        end
    end
end


text:
function ClassificationDPL(TtData::Matrix{Float64}, DictMat::Array{Any,1}, EncoderMat::Matrix{Float64}, DictSize::Int64)

    # Projective representation coefficients estimation
    ClassNum::Int64              = length(DictMat)
    PredictCoef::Matrix{Float64} = EncoderMat * TtData

    Error::Matrix{Float64} = zeros(ClassNum, size(TtData, 2));
    # Class-specific reconstruction error calculation
    for i=1:ClassNum
        @inbounds reconstructedTtData::Matrix{Float64} = DictMat[i] * PredictCoef[(i-1)*DictSize+1:i*DictSize, :]
        reconstructedTtData -= TtData
        @inbounds Error[i,:] = sum(abs2(reconstructedTtData), 1)
    end
    Distance::Matrix{Float64}, PredictInd::Matrix{Int64} = findmin(Error, 1)
    PredictLabel = [ind2sub(size(Error), PredictInd[i])[1] for i = 1:size(PredictInd, 2)]

    PredictLabel, Error
end


text:
@testset "VehicleState" begin

    # Test Base
    x = VehicleState()
    xc = copy(x)
    @test x == xc
    xc.x = 10
    @test x != xc

    x1 = VehicleState(1.0,2.0,3.0,4.0)
    x2 = VehicleState(2.0,3.0,4.0,5.0)
    x3 = x1 + x2
    @test x3.x == 3.0
    @test x3.y == 5.0
    @test x3.θ == 7.0
    @test x3.v == 9.0

    x4 = x2*0.5
    @test x4.x == 1.0
    @test x4.y == 1.5
    @test x4.θ == 2.0
    @test x4.v == 2.5

    x4 = 0.5*x2
    @test x4.x == 1.0
    @test x4.y == 1.5
    @test x4.θ == 2.0
    @test x4.v == 2.5

    # Test standardize
    # DoubleIntegrator
    model = Algames.DoubleIntegratorGame(p=3)
    xi = @SVector [1.0, 2.0, sqrt(3)/2, 0.5]
    vsi = standardize(model, xi)

    @test vsi.x == xi[1]
    @test vsi.y == xi[2]
    @test abs(vsi.θ - pi/6) < 1e-10
    @test abs(vsi.v - 1) < 1e-10

    x = @SVector [1.0, 1.2, 1.4,
                  2.0, 2.2, 2.4,
                  sqrt(3)/2, sqrt(3)/2, sqrt(3)/2,
                  0.5, 0.5, 0.5]
    vs = standardize(model, x)

    @test vs[1].x == 1.0
    @test vs[1].y == 2.0
    @test abs(vs[1].θ - pi/6) < 1e-10
    @test abs(vs[1].v - 1) < 1e-10

    @test vs[2].x == 1.2
    @test vs[2].y == 2.2
    @test abs(vs[2].θ - pi/6) < 1e-10
    @test abs(vs[2].v - 1) < 1e-10

    @test vs[3].x == 1.4
    @test vs[3].y == 2.4
    @test abs(vs[3].θ - pi/6) < 1e-10
    @test abs(vs[3].v - 1) < 1e-10

    # Unicycle
    model = UnicycleGame(p=3)
    xi = @SVector [1.0, 2.0, sqrt(3)/2, 0.5]
    vsi = standardize(model, xi)

    @test vsi.x == xi[1]
    @test vsi.y == xi[2]
    @test vsi.θ == xi[3]
    @test vsi.v == xi[4]

    # Bicycle
    model = BicycleGame(p=3)
    xi = @SVector [1.0, 2.0, sqrt(3)/2, 0.5]
    vsi = standardize(model, xi)

    @test vsi.x == xi[1]
    @test vsi.y == xi[2]
    @test vsi.θ == xi[4]
    @test vsi.v == xi[3]

    # Test specialize
    T = Float64
    p = 3
    model = DoubleIntegratorGame(p=p)
    x = rand(SVector{model.n,T})
    vs = standardize(model, x)
    x_ = specialize(model, vs)
    @test norm(x - x_, 1) < 1e-10

    p = 3
    model = UnicycleGame(p=p)
    x = rand(SVector{model.n,T})
    vs = standardize(model, x)
    x_ = specialize(model, vs)
    @test norm(x - x_, 1) < 1e-10

    p = 3
    model = BicycleGame(p=p)
    x = rand(SVector{model.n,T})
    vs = standardize(model, x)
    x_ = specialize(model, vs)
    @test norm(x - x_, 1) < 1e-10

end


text:
function fib10()
#
#  test polynomial suggested by Goedecker
#
    n = 10
    p = fib(n)
    z = [-.8990310941957183    1;
         -.7399333995096503-.5167965365015562*im  1;
         -.7399333995096503+.5167965365015562*im  1;
         -.3130259587265545-.8583547482505247*im  1;
         -.3130259587265545+.8583547482505247*im  1;
         .2462366020696891-.9013280059509374*im  1;
         .2462366020696891+.9013280059509374*im  1;
         .7567289869093243-.6038877405652167*im  1;
         .7567289869093243+.6038877405652167*im  1;
         1.999018632710101  1; ]
    p, PolyZeros(z)
end


text:
class A{
  method() {
    print "A method";
  }
}

class B < A {
  method() {
    print "B method";
  }

  test() {
    super.method();
  }
}

class C < B {}

C().test();

class Doughnut {
  cook() {
    print "Fry until golden brown.";
  }
}

class BostonCream < Doughnut {
  cook() {
    super.cook();
    print "Pipe full of custard and coat with chocolate.";
  }
}

BostonCream().cook();


text:
module MatrixNetworks
	using LinearAlgebra
	using DelimitedFiles

	export random_adjmat, symmetrize, neighbours, degree, load_matrix_network

	"""
		symmetrize(M::Matrix) -> SymM::Symmetric

	Convert rectangular (`n_L` × `n_T`) matrix into a symmetric (`n_L` + `n_T`) order
	matrix.

	# Examples
	```julia-repl
	julia> A = [ 1 0 1 ; 0 0 1 ]
	2×3 Array{Int64,2}:
	1  0  1
	0  0  1
	julia> symmetrize(A)
	5×5 Array{Int64,2}:
	0  0  1  0  1
	0  0  0  0  1
	1  0  0  0  0
	0  0  0  0  0
	1  1  0  0  0
	```

	# Arguments
	- `M::Matrix`: rectangular (`n_L` × `n_T`) matrix.

	"""
	function symmetrize(M::Matrix)::Symmetric
		# Get dimensions of non-symmetric matrix
		R, C = size(M)

		# Create empty matrices to fill spaces between our matrix
		Rzeros = zeros(Int64, R, R)
		Czeros = zeros(Int64, C, C)

		# Create symmetric matrix
		M	 = [ Rzeros M ; M' Czeros ]
		SymM = Symmetric(M)

		return SymM
	end

	"""
		neighbours(M::Matrix, node::Int, orientation::String = "row") -> neighbours_array::Array{Int}

	Get neighbours of node in adjacency matrix.

	# Examples
	```julia-repl
	julia> A = [ 0 0 1 0 1 ; 0 0 0 0 1 ; 1 0 0 0 0 ; 0 0 0 0 0 ; 1 1 0 0 0 ]
	5×5 Array{Int64,2}:
	0  0  1  0  1
	0  0  0  0  1
	1  0  0  0  0
	0  0  0  0  0
	1  1  0  0  0
	julia> neighbours(A, 1)
	2-element Array{Int64,1}:
	3
	5
	```

	```julia-repl
	julia> B = [ 0 1 0 0 ; 1 0 1 0 ]
	2×4 Array{Int64,2}:
	0  1  0  0
	1  0  1  0
	julia> neighbours(B, 2, "row")
	2-element Array{Int64,1}:
	1
	3
	julia> neighbours(B, 2, "column")
	1-element Array{Int64,1}:
	2
	```

	# Arguments
	- `M::Matrix`: adjacency matrix.
	- `node::Int`: idx of node.
	- `orientation::String="row"` : orientation of search (column or row-wise).

	"""
	function neighbours(M::Matrix, node::Int, orientation::String = "row")::Array{Int}
		if orientation == "row"
			node_adjacency	 = M[node,:]
		elseif orientation == "column"
			node_adjacency	 = M[:, node]
		else
			error("Orientation value must be either \"row\" of \"column\"")
		end

		neighbours_array = [idx for (idx,neighbour) in enumerate(node_adjacency)
							if neighbour != 0]

		return neighbours_array
	end

	"""
		degree(M::Matrix, node::Int, orientation::String = "row") -> node_degree::Int

	Get degree of node in symmetric adjacency matrix.

	# Examples
	```julia-repl
	julia> A = [ 0 0 1 0 1 ; 0 0 0 0 1 ; 1 0 0 0 0 ; 0 0 0 0 0 ; 1 1 0 0 0 ]
	5×5 Array{Int64,2}:
	0  0  1  0  1
	0  0  0  0  1
	1  0  0  0  0
	0  0  0  0  0
	1  1  0  0  0
	julia> degree(A, 1)
	2
	```

	```julia-repl
	julia> B = [ 0 1 0 0 ; 1 0 1 0 ]
	2×4 Array{Int64,2}:
	0  1  0  0
	1  0  1  0
	julia> degree(B, 2, "row")
	2
	julia> degree(B, 2, "column")
	1
	```

	# Arguments
	- `M::Matrix`: adjacency matrix.
	- `node::Int`: idx of node.
	- `orientation::String="row"` : orientation of search (column or row-wise).

	"""
	function degree(M::Matrix, node::Int, orientation::String = "row")::Int
		if orientation == "row"
			node_neighbours = neighbours(M, node, "row")
		elseif orientation == "column"
			node_neighbours = neighbours(M, node, "column")
		else
			error("Orientation value must be either \"row\" of \"column\"")
		end

		node_degree = length(node_neighbours)

		return node_degree
	end

	"""
		load_matrix_network(filepath::String, separator::Char = ' ') -> A::Matrix{Int64}
	"""
	function load_matrix_network(filepath::String, separator::Char = ' ')::Matrix{Int64}
		A = readdlm(filepath, separator, Int64, '\n')
		return A
	end

end


text:
# -------------------------------------------------------------------------------------
# _SEPS
const _SEPS = Dict{Symbol, Char}()

function _filesep()
    dummy = joinpath("a", "a")
    reg = Regex("^a(?<sep>.*)a\$")
    sep = match(reg, dummy)["sep"]
    (length(sep) != 1) && error("Unsupported file separator")
    first(sep)
end

function _set_default_SEPS!()
    empty!(_SEPS)
    _SEPS[:EXT_SEP] = '.'
    _SEPS[:PAIR_SEP] = '='
    _SEPS[:ELEMT_SEP] = ' '
    _SEPS[:PARAMS_LSEP] = '['
    _SEPS[:PARAMS_RSEP] = ']'
    _SEPS[:SEP_SUBST] = '_'
    _SEPS[:FILE_SEP] = _filesep()
    _check__SEPS()
    return _SEPS
end

function _check__SEPS()
    isempty(_SEPS) && error("_SEPS is empty")
    !allunique(values(_SEPS)) && error("_SEPS are not unique. Current _SEPS ", _SEPS)
    return _SEPS
end

_hex_escaped_seps(ks = keys(_SEPS)) = Dict{Symbol, String}(k => hex_escape(_SEPS[k]) for k in ks)

# -------------------------------------------------------------------------------------
# RESERVED SEPS
const _RESERVED_SEPS_KEYS = [:ELEMT_SEP, :PAIR_SEP, :PARAMS_LSEP, :PARAMS_RSEP, :FILE_SEP]
_reserved_seps() = [_SEPS[S] for S in _RESERVED_SEPS_KEYS]

# -------------------------------------------------------------------------------------
# INPUT TYPES
const _INPUT_KEY_TYPES = [Symbol, String]
const _INPUT_VAL_TYPES = [Symbol, Bool, AbstractString, Integer, AbstractFloat]
const _INPUT_PAIRS_TYPES = [Pair, Dict, NamedTuple]

_iskeyT(k) = any(isa.([k], _INPUT_KEY_TYPES))
_isvalT(v) = any(isa.([v], _INPUT_VAL_TYPES))
_ispairT(p) = any(isa.([p], _INPUT_PAIRS_TYPES))

function _checker(v, isfun, name, types)
    !isfun(v) &&
        error("Invalid input ", name, " type, got ", v , "::",  typeof(v),
            ", expected a type from [", join(types, ", "), "]"
        )
    return v
end

_check_keyT(k) = _checker(k, _iskeyT, "key", _INPUT_KEY_TYPES)
_check_valT(v) = _checker(v, _isvalT, "value", _INPUT_VAL_TYPES)
_check_pairT(p) = _checker(p, _ispairT, "key:value structure", _INPUT_PAIRS_TYPES)

# -------------------------------------------------------------------------------------
# String validity
const INVALID_CHARS = Char[]

function _check_str(str::String)

    reserved_seps = _reserved_seps()
    any(contains.(str, reserved_seps)) &&
        error("Separator detected at '$str'. Reserved seps $reserved_seps")
    any(contains.(str, INVALID_CHARS)) &&
        error("Invalid char detected '$str'. INVALID_CHARS $INVALID_CHARS")

    return str
end

# -------------------------------------------------------------------------------------
const _OUTPUT_VAL_TYPES = [Int, Float64, Bool, String] # String must be the last, it is the fallback


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

mutable struct InfoField <: SwaggerModel
    project::Any # spec type: Union{ Nothing, String } # spec name: project

    function InfoField(;project=nothing)
        o = new()
        validate_property(InfoField, Symbol("project"), project)
        setfield!(o, Symbol("project"), project)
        o
    end
end # type InfoField

const _property_map_InfoField = Dict{Symbol,Symbol}(Symbol("project")=>Symbol("project"))
const _property_types_InfoField = Dict{Symbol,String}(Symbol("project")=>"String")
Base.propertynames(::Type{ InfoField }) = collect(keys(_property_map_InfoField))
Swagger.property_type(::Type{ InfoField }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_InfoField[name]))}
Swagger.field_name(::Type{ InfoField }, property_name::Symbol) =  _property_map_InfoField[property_name]

function check_required(o::InfoField)
    true
end

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


text:
using OpenQuantumTools, LaTeXStrings
using OrdinaryDiffEq, Plots, StatsBase

# calculate the mean and standard deviation of the mean estimator from a sample
function mean_std(sample)
    m, v = mean_and_std(sample, corrected=true)
    m, v/sqrt(length(sample))
end

# All values calculated below are in angular frequency units
num = 10
bvec = 0.2 * ones(num)
γvec = log_uniform(0.01, 1, num)
fluctuator_ensemble = EnsembleFluctuator(bvec, γvec);

plot(fluctuator_ensemble, :spectrum, 2*π*range(0.01, 4, length=100), xscale=:log10, yscale=:log10, linewidth=2, label="")
xlabel!(L"\omega")
ylabel!(L"S(\omega)")

H = DenseHamiltonian([(s)->1], [σz], unit=:ħ)
coupling = ConstantCouplings([0.5*σz], unit=:ħ)
u0 = PauliVec[1][1]
annealing = Annealing(H, u0, coupling = coupling, bath=fluctuator_ensemble)
tf = 10
# create object for parallel simulation
prob = build_ensembles(annealing, tf, :stochastic)
# we run each trajectory serially for this example
sol = solve(prob, Tsit5(), EnsembleSerial(), trajectories=1000, reltol=1e-6, saveat=range(0,tf,length=100))

t_axis = range(0,tf,length=100)
es = []
err = []
for (i, s) in enumerate(t_axis)
    sample = [abs2(u0'*so[i]) for so in sol]
    pop, pop_std = mean_std(sample)
    push!(es, pop)
    push!(err, 2*pop_std)
end

plot(t_axis, es, ribbon=err, linewidth=2, label="")
xlabel!("ns")
ylabel!("<+>")

cb = InstPulseCallback([0.25, 0.5, 0.75] .* tf, (c, x) -> c .= σx * c)
sol = solve(prob, Tsit5(), EnsembleSerial(), trajectories=1000, reltol=1e-6, saveat=range(0,tf,length=100), callback=cb)

s_axis = range(0,1,length=100)
es = []
err = []
for (i, s) in enumerate(s_axis)
    sample = [abs2(u0'*so[i]) for so in sol]
    pop, pop_std = mean_std(sample)
    push!(es, pop)
    push!(err, 2*pop_std)
end

plot(tf*s_axis, es, ribbon=err, linewidth=2, label="")
xlabel!("ns")
ylabel!("<+>")

using HOQSTTutorials
HOQSTTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
"""
Module with general utilities and missing standard library features
that could be useful in any Julia project
"""
module Utils

include("misc.jl")
include("information.jl")
include("sample.jl")

end #module


text:
# To run:
# include("DI_design.jl")

# DESCRIPTION
# Trims and linearizes the rotorcraft model at incremental flight speeds, and
# generates the gains of the dynamic inversion controller. Always run after
# modifying the aircraft parameters in the AIRCRAFTNAME_constants.jl file.
#
# ------------------------------------------------------------------------------

# include packages
#using DifferentialEquations
using LinearAlgebra
using SparseArrays
using ControlSystems
using MAT
using Interpolations
using ApproxFun
using Printf
using Plots
using JLD
using LaTeXStrings

# include functions
include("atan2.jl")
include("linearize.jl")
include("trimmer.jl")
include("SimpleHel.jl")
#include("rk4.jl")
#include("finp.jl")
# load H-60 constants
include("H60_constants.jl")
# vector of trim speeds [kts]
VXTRIM_vec=[0.01,20.0,40.0,60.0,80.0,100.0,120.0,140.0,160.0]
# initialize DI matrices
CATAB = zeros(3,3,length(VXTRIM_vec))
CBinvTAB = zeros(3,3,length(VXTRIM_vec))
x0_mat = zeros(NSTATES,length(VXTRIM_vec))
u0_mat = zeros(NCTRLS,length(VXTRIM_vec))
# trim at different speeds
for iv=1:length(VXTRIM_vec)
    # load H60 constants
    include("H60_constants.jl")
    # forward, lateral, and velocities in trim [ft/s]
    VXTRIM=VXTRIM_vec[iv]*1.688
    VYTRIM=0.0
    VZTRIM=0.0
    # turn rate in trim [rad/s]
    PSIDTRIM=0.0
    # initial heading [rad]
    PSITRIM = 0.0*pi/180
    # initial position [ft]
    ALTTRIM = 0.0
    XNTRIM = 0.0
    YETRIM = 0.0
    # initialize state vectors and initial guess
    x0=zeros(NSTATES)
    x0[1]=VXTRIM
    x0[2]=VYTRIM
    x0[3]=VZTRIM
    # option to trim in zero sideselip or zero bank angle
    if (VXTRIM>=60*1.688)
        if PSIDTRIM==0
            # Trim heading for zero bank angle above 60 kts
            x0[7]=0.0;
            TRIMVARS[7:8]=[8; 9]
        else
            x0[9]=atan2(VYTRIM,VXTRIM)
            TRIMVARS[7:8]=[7; 8]
            # initial guess for phi (using turn coord)
            x0[7]=atan(PSIDTRIM*sqrt(VXTRIM^2+VYTRIM^2+VZTRIM^2)/G)
            # initial guess for r
            x0[6]=G/sqrt(VXTRIM^2+VYTRIM^2+VZTRIM^2)*sin(x0[7])
            # initial guess for q
            x0[5]=x0[6]*tan(x0[7])
        end
    else
        # otherwise set heading to flight path
        x0[9]=atan2(VYTRIM,VXTRIM)
        TRIMVARS[7:8]=[7; 8]
    end
    # initial guess for controls [deg]
    u0=[50.0; 50.0; 50.0; 50.0]
    # target values for state derivatives
    xdot_targ=zeros(12)
    xdot_targ[9]=PSIDTRIM
    xdot_targ[10]=VXTRIM
    xdot_targ[11]=VYTRIM
    xdot_targ[12]=VZTRIM
    # time [s]
    t=0.0
    # trim aircraft
    x0, u0 = trimmer(SimpleHel!,x0,u0,xdot_targ)
    # linearized aircraft dynamics
    A, B, C, D = linearize(SimpleHel!,x0,u0)
    # store trim state and controls
    x0_mat[:,iv] = x0
    u0_mat[:,iv] = u0
    # dynamic inverse matrices
    C_DI=[1 0 0
          0 1 0
          0 0 1]
    CATAB[:,:,iv]=C_DI*A[4:6,4:6]
    CBinvTAB[:,:,iv]=inv(C_DI*B[4:6,[1,2,4]])
end

# command filters time constants
tau_p=1/3.5
tau_q=1/4.5
tau_r=0.5
# washout filter for controls
tau_u=5.0;
# disturbance rejection frequency and damping and integrator poles
wnroll_d=3.5
dmproll_d=1.0
wnpitch_d=4.5
dmppitch_d=1.0
wnyaw_d=2.;
dmpyaw_d=1.;
# feedback gains
kp_roll=2*wnroll_d*dmproll_d
ki_roll=wnroll_d^2
kp_pitch=2*wnpitch_d*dmppitch_d
ki_pitch=wnpitch_d^2
kp_yaw=2*dmpyaw_d*wnyaw_d
ki_yaw=wnyaw_d^2;
# pilot gains
k_lat=pi/100
k_lon=pi/100
k_ped=0.644
# save DI matrices and gains into a file
save("DI.jld","CATAB",CATAB,"CBinvTAB",CBinvTAB,"kp_roll",kp_roll,"ki_roll",
    ki_roll,"kp_pitch",kp_pitch,"ki_pitch",ki_pitch,"kp_yaw",kp_yaw,"ki_yaw",
    ki_yaw,"k_lat",k_lat,"k_lon",k_lon,"k_ped",k_ped,"tau_p",tau_p,
    "tau_q",tau_q,"tau_r",tau_r,"tau_u",tau_u)

# store trim state and control input vectors
save("trim_J-SimpleHel.jld","x0_mat",x0_mat,"u0_mat",u0_mat,"VXTRIM_vec",VXTRIM_vec)

# plot trim attitudes
#gr(size=(1000,1000))
gr(size=(800,600))
p1=plot(VXTRIM_vec,x0_mat[7,:]*180/pi,label="J-GenHel",legend=:bottomright, legendfontsize=12,line=(3, :blue, :solid), marker = (:circ, :blue, 8), markerstrokecolor = :blue)
xaxis!((0, 160), 0:20:160)
yaxis!("φ [deg]",yguidefontsize=14)
p2=plot(VXTRIM_vec,x0_mat[8,:]*180/pi,label="",line=(3, :blue, :solid), marker = (:circ, :blue, 8), markerstrokecolor = :blue)
xaxis!((0, 160), 0:20:160)
yaxis!("θ [deg]",yguidefontsize=14)
p3=plot(VXTRIM_vec,x0_mat[9,:]*180/pi,label="",line=(3, :blue, :solid), marker = (:circ, :blue, 8), markerstrokecolor = :blue)
xaxis!((0, 160), 0:20:160)
xaxis!("Total speed, V [kts]",xguidefontsize=14)
yaxis!("ψ [deg]",yguidefontsize=14)
p4=plot(p1, p2, p3, layout=(3,1))
display(p4)
#savefig(".\\Plots\\trim.svg")

# plot trim controls
gr(size=(800,600))
p5=plot(VXTRIM_vec,u0_mat[1,:],label="J-SimpleHel",legend=:bottomright, legendfontsize=12,line=(3, :blue, :solid), marker = (:circ, :blue, 8), markerstrokecolor = :blue)
xaxis!((0, 160), 0:20:160)
yaxis!(L"\delta_\mathrm{lat} [\%]",yguidefontsize=14)
p6=plot(VXTRIM_vec,u0_mat[2,:],label="",line=(3, :blue, :solid), marker = (:circ, :blue, 8), markerstrokecolor = :blue)
xaxis!((0, 160), 0:20:160)
yaxis!(L"\delta_\mathrm{lon} [\%]",yguidefontsize=14)
p7=plot(VXTRIM_vec,u0_mat[3,:],label="",line=(3, :blue, :solid), marker = (:circ, :blue, 8), markerstrokecolor = :blue)
xaxis!((0, 160), 0:20:160)
yaxis!(L"\delta_\mathrm{col} [\%]",yguidefontsize=14)
p8=plot(VXTRIM_vec,u0_mat[4,:],label="",line=(3, :blue, :solid), marker = (:circ, :blue, 8), markerstrokecolor = :blue)
xaxis!((0, 160), 0:20:160)
xaxis!("Total speed, V [kts]",xguidefontsize=14)
yaxis!(L"\delta_\mathrm{ped} [\%]",yguidefontsize=14)
p9=plot(p5, p6, p7, p8, layout=(4,1))
display(p9)
#savefig(".\\Plots\\trim_ctrls.svg")

max_stars_repo_name: umbertosaetti/J-SimpleHel

max_issues_repo_name: umbertosaetti/J-SimpleHel

max_forks_repo_name: umbertosaetti/J-SimpleHel


text:
todisplay = VERSION <= VersionNumber("1.2") ? true : false

using LinearAlgebra
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using CAGD
if todisplay
    using ViewerGL
    GL = ViewerGL
	include("views.jl")
end

# Rod Generation

# npts = 4
npts = 16
V,_ = Lar.rod()([npts, 1])
EV = [
    [[2*i, (2*i+1)%2npts+1] for i = 1 : npts];             # horizontal upper edges
    [[2*i-1, (2*i)%2npts+1] for i = 1 : npts];           # horizontal lower edges
    [[2*i+1, 2*i+2] for i = 0 : npts-1];                # vertical edges
]
FE = [
    [i for i = 1 : npts],                                       # upper face
    [npts + i for i = 1 : npts],                                # lower face
    [[i, npts+i, 2*npts+i, 2*npts+i+1] for i = 1:npts-1]...,    # vertical faces
    [npts, 2*npts, 3*npts, 2*npts+1]
]
FV = Lar.cop2lar(Lar.lar2cop(FE) * Lar.lar2cop(EV))
CV = [collect(1 : size(V,2))]

if todisplay
    GL.VIEW([
        GL.GLGrid(V,EV,GL.COLORS[1],1)
        GL.GLFrame
    ]);
end

# Three Rods Generation

cyl = Lar.Struct([
    Lar.r(0,0,0),
    #Lar.s(0.6, 0.6, 2.0),
    Lar.s(1.0, 1.0, 1.5),
    Lar.t(0,0,-1.5),
    (V,CV,FV,EV)
])
cyl = Lar.struct2lar(cyl)
tris = Lar.Struct([
    cyl,
    Lar.Struct([Lar.r(pi/2,0,0), cyl ]) ,
    Lar.Struct([Lar.r(0,pi/2,0), cyl ])
])
V, CV, FV, EV = Lar.struct2lar(tris)
cyls = (V, CV, FV, EV)

# Cube Generation

V, (_, EV, FV, CV) = Lar.cuboid([1,1,1],true,[-1,-1,-1])
cube = (V, CV, FV, EV)

# Sphere Generation

# V, EV, FV = catmullclark(cube[1], cube[3], cube[2], 4)
# sphere = (V, FV, EV)

V, FV = Lar.sphere(1.0)([15,25])
# V, FV = Lar.sphere(1.0)([5,5])
FV = sort(sort.(FV))
function getFaceEdges(fV::Array{Int,1})
    return [fV[1], fV[2]], [fV[1], fV[3]], [fV[2], fV[3]]
end
EVs = unique([(map(f -> getFaceEdges(f), FV)...)...])
CV = [collect(1:size(V, 2))]
sphere = (V, CV, FV, EVs)

# Object Creation

carry = Lar.Struct([
    cyls,
    Lar.s(1.3,1.3,1.3),
    cube,
    Lar.s(1.4,1.4,1.4),
    #Lar.s(1.5,1.5,1.5)
    sphere
])
V,CV,FV,EV = Lar.struct2lar(carry)

if todisplay
    GL.VIEW([
        GL.GLGrid(V,EV,GL.COLORS[1],1)
        GL.GLFrame
    ]);
end

# Model Generation

model = CAGD.Model(V)
cFE = convert(Lar.ChainOp, Lar.coboundary_1(model.G, FV, EV))
CAGD.addModelCells!(model, 1, EV, signed = true)
CAGD.addModelCells!(model, 2, cFE)
# Since CF is only used in boolean evaluation, the sign is not needed:
cFV = Lar.lar2cop(FV)
cCV = Lar.lar2cop(CV)
cCF = convert(Lar.ChainOp, ((cCV*cFV').>0))
CAGD.addModelCells!(model, 3, cCF)

atol = 1e-9;
if todisplay  displayModel(model)  end

split_model = CAGD.facesplitting(model, atol = atol)

if todisplay  displayModel(split_model)  end

congr_model = CAGD.mergeModelVertices(split_model, err=atol, signed_merge=true)
#split_model = CAGD.facesplitting(congr_model, atol = atol)
#congr_model = CAGD.mergeModelVertices(split_model, err=atol, signed_merge=true)

if todisplay  displayModel(congr_model)  end

gift_model = deepcopy(congr_model);
FC, bicon_comps = CAGD.tgw(congr_model, 3)
CAGD.addModelCells!(gift_model, 3, convert(Lar.ChainOp, FC'))

if todisplay  viewExplode(gift_model)  end

##==============================================================================
##  Boolean Decomposition
##==============================================================================

arranged_model, boolean_matrix = CAGD.bool3(model)

bXRod = boolean_matrix[:, 2]
bYRod = boolean_matrix[:, 3]
bZRod = boolean_matrix[:, 4]
bCube = boolean_matrix[:, 5]
bSphe = boolean_matrix[:, 6]


text:
using githubactions4julia
using Test

@testset "githubactions4julia.jl 基本テスト" begin
    # Write your tests here.
    @test githubactions4juliasample()
end