text:
# Analyze Monsoon rainfail data from
# @Article{sinha2015trends,
#   Title                    = {Trends and oscillations in the Indian summer monsoon rainfall over the last two millennia},
#   Author                   = {Sinha, A. and Kathayat, G. and Cheng, H. and others},
#   Journal                  = {Nature communications},
#   Year                     = {2015},
#   Volume                   = {6},
#   Publisher                = {Nature Publishing Group},
# }
using EstimateSDE
using Plots
using DataFrames
using Pipe:@pipe
pgfplotsx()
using Serialization
using Parameters
using Printf

## Load results rainfall-estimate
results_dir = EstimateSDE.resultsdir("rainfall-data")
f = joinpath(results_dir,"rainfall-sde-estimate.jld")
R = deserialize(f); @info "Results loaded to: $f"

## Settings
nspec = 500 # Number of points in the spectrum
ssize = 300,200

@unpack models,settings = R
@unpack Tkld,Tsub = settings

colors = [:blue,:red,:green,:orange,:yellow]

## Report
fmax = 1/2Tkld
f = (0:nspec-1)/nspec*1/2Tkld
plt_psd = plot(xlabel="f [1/year]",ylabel="psd",size=ssize,fg_legend=nothing,bg_legend=nothing)
for (i,m) in enumerate(models[1:end-1])
  h = powerspectrum(m.proc,f)
  kwargs = if m.label=="original"
    (color=:black,linestyle=:dash,label="all data")
  else
    (linealpha=0.8,label=nothing,color=colors[i-1])
  end
  plot!(f,h,label=m.label,yscale=:log10;kwargs...)
  println(m.label)
  describe_model(m.proc)
end
fmax_sub = 1/2Tsub
plot!([fmax_sub],linestyle=:dash,linecolor=:black,seriestype=:vline,label=nothing,
    legend=:topright)
# display(plt_psd) # fails in vscode with linealpha problem
@info "PSD plot not displayed, because it results in error in VS Code"

## trace plot
t_kld = 0:Tkld:100
proc = models[1].proc
plt_order = plot(xlabel="order p",ylabel="-L",fg_legend=nothing,bg_legend=nothing,size=ssize,
    legend=:topright)
for (i,m) in enumerate(models)
  if m.label=="original" continue end
  trace = m.trace
  data = m.data
  R,Href = tracetable(data,trace,proc,t_kld)

  # Maximum Likelihood per order
  mlo = @pipe(
    groupby(R,:p) |>
    combine(_) do df
      df[argmax(df.L),:]
    end |>
    sort(_,:p)
  )
  mlo_cols = [:p,:Hhat,:D,:init_tor_type,:stage_str]
  display(mlo[:,mlo_cols])
  # Overall maximum likelihood
  Imle = R[argmax(R.L),:]

  ## Plot trace
  Hhat_min = minimum(mlo.Hhat)
  Hhat = map(x->isinf(x) ? missing : x,mlo.Hhat)
  plot!(mlo.p,Hhat.-Hhat_min,ms=2,msc=:auto,markershape=:rect,label=nothing,
    color=colors[i-1],linealpha=0.8)
  plot!([Imle.p],[Imle.Hhat-Hhat_min],lc=nothing,mc=:yellow,markershape=:diamond,ma=0.5,ms=8,
      label=i==2 ? "ML" : nothing)
end
# @info "Plot fit vs order not shown"
plot!(ylim=(0,8))
# display(plt_order)

plt = plot(plt_psd,plt_order,layout=(2,1),size=(ssize[1],2ssize[2]))

## Save figure
for ext in ["pdf","tex"]
  f = joinpath(results_dir,"rainfall-results.$ext")
  savefig(plt,f)
  @info "Figure saved to $f"
end


text:
# translator for the APS2
module APS2

import Base: convert

using HDF5

import QGL

import ..config

const DAC_CLOCK = 1.2e9
const FPGA_CLOCK = 300e6
const ADDRESS_UNIT = 4  #everything is done in units of 4 timesteps
const MIN_ENTRY_LENGTH = 8
const MAX_WAVEFORM_PTS = 2^28  #maximum size of waveform memory
const WAVEFORM_CACHE_SIZE = 2^17
const MAX_WAVEFORM_VALUE = 2^13 - 1  #maximum waveform value i.e. 14bit DAC
const MAX_NUM_INSTRUCTIONS = 2^26
const MAX_REPEAT_COUNT = 2^16 - 1
const MAX_MARKER_COUNT = 2^32 - 1

# instruction encodings
const WFM = 0x00
const MARKER = 0x01
const WAIT = 0x02
const LOAD_REPEAT = 0x03
const DEC_REPEAT = 0x04
const CMP = 0x05
const GOTO = 0x06
const CALL = 0x07
const RETURN = 0x08
const SYNC = 0x09
const MODULATOR = 0x0a
const LOAD_CMP = 0x0b
const PREFETCH = 0x0c

const APS2Instruction = UInt64

immutable Waveform
	address::UInt32
	count::UInt32
	isTA::Bool
	write_flag::Bool
	instruction::UInt64
end

# WFM/MARKER op codes
const PLAY_WFM      = 0x0
const WAIT_TRIG     = 0x1
const WAIT_SYNC     = 0x2
const WFM_PREFETCH  = 0x3
const WFM_OP_OFFSET = 46
const TA_PAIR_BIT   = 45
const WFM_CT_OFFSET = 24

function Waveform(address, count, isTA, write_flag)
	ct = UInt64(count ÷ ADDRESS_UNIT - 1) & 0x000f_ffff # 20 bit count
	addr = UInt64(address ÷ ADDRESS_UNIT) & 0x00ff_ffff # 24 bit address
	header = UInt64( (WFM << 4) | (0x3 << 2) | (write_flag & 0x1) )
	payload = (UInt64(PLAY_WFM) << WFM_OP_OFFSET) | (UInt64(isTA) << TA_PAIR_BIT) | (ct << WFM_CT_OFFSET) | addr
	instr = (header << 56) | payload
	Waveform(addr, ct, isTA, write_flag, instr)
end

immutable Marker
	engine_select::UInt8
	state::Bool
	count::UInt32
	transition_word::UInt8
	write_flag::Bool
	instruction::UInt64
end

function Marker(marker_select, count, state, write_flag)
	count = UInt64(count)
	quad_count =  UInt64(count) ÷ ADDRESS_UNIT & UInt64(0x0fff_ffff) # 28 bit count
	count_rem = count % ADDRESS_UNIT
	if state
		transition_words = [0b1111; 0b0111; 0b0011; 0b0001]
		transition_word = transition_words[count_rem+1]
	else
		transition_words = [0b0000; 0b1000; 0b1100; 0b1110]
		transition_word = transition_words[count_rem+1]
	end
	header = (MARKER << 4) | (((marker_select-1) & 0x3) << 2) | (write_flag & 0x1)
	payload = (UInt64(PLAY_WFM) << WFM_OP_OFFSET) | (UInt64(transition_word) << 33) | (UInt64(state) << 32) | quad_count
	instr = (UInt64(header) << 56) | payload
	Marker(UInt8(marker_select), state, quad_count, transition_word, write_flag, instr)
end

immutable ControlFlow
	instruction::UInt64
end

# modulation instructions
@enum MODULATION_OP_CODE MODULATE=0x00 RESET_PHASE=0x02 SET_FREQ=0x06 SET_PHASE=0x0a UPDATE_FRAME=0x0e
const MODULATOR_OP_OFFSET = 44
const NCO_SELECT_OP_OFFSET = 40
function modulation_instr(op::MODULATION_OP_CODE, nco_select, payload=0)
	return UInt64(MODULATOR) << 60 | UInt64(0x1) << 56 |
		UInt64(op) << MODULATOR_OP_OFFSET | UInt64(nco_select) << 40 |
		reinterpret(UInt32, Int32(payload))
end

"""
Serialize a pulse sequence to a HDF5 file
"""
function write_sequence_file(filename, seqs, pulses, channel_map)

	# check whether there is any analog data
	markers_only = !(:ch12 in keys(channel_map))
	# translate pulses to waveform and/or markers
	instr_lib = Dict{QGL.Pulse, Union{Waveform,Marker}}()
	chan_freqs = Dict{QGL.Channel, Float64}()
	wfs = Vector{Vector{Complex{Int16}}}()
	if !markers_only
		for ch in channel_map[:ch12]
			chan_freqs[ch] = ch.frequency
			create_wf_instrs!(instr_lib, wfs, pulses[ch])
		end
	end
	for (ct, marker_chan) = enumerate([:m1, :m2, :m3, :m4])
		if marker_chan in keys(channel_map)
			create_marker_instrs!(instr_lib, pulses[channel_map[marker_chan][1]], ct)
		end
	end

	# create instructions
	instrs = create_instrs(seqs, instr_lib, collect(values(channel_map)), chan_freqs)

	write_to_file(filename, instrs, wfs)
end

const USE_PHASE_OFFSET_INSTRUCTION = false
const USE_PULSE_FREQUENCY_INSTRUCTION = false

function create_wf_instrs!(instr_lib, wfs, pulses)
	# TODO: better handle Id so we don't generate useless long wfs and have repeated TAZ offsets
	idx = 0
	for p in pulses
		wf = p.amp * QGL.waveform(p, DAC_CLOCK)
		if !USE_PHASE_OFFSET_INSTRUCTION
			wf *= exp(1im * 2π * p.phase)
		end
		if !USE_PULSE_FREQUENCY_INSTRUCTION && p.frequency != 0
			# bake the pulse frequency into the waveform
			wf .*= exp.(-1im * 2π * p.frequency * (1/DAC_CLOCK) * (1:length(wf)) )
		end
		# reduce to Int16 with maximum for 14 bit DAC
		wf = round.(Int16, MAX_WAVEFORM_VALUE*real(wf)) + 1im*round.(Int16, MAX_WAVEFORM_VALUE*imag(wf))

		isTA = all(wf .== wf[1])
		instr_lib[p] = Waveform(idx, length(wf), isTA, true)
		if isTA
			idx += ADDRESS_UNIT
			push!(wfs, wf[1:ADDRESS_UNIT])
		else
			idx += length(wf)
			push!(wfs, wf)
		end
	end

	return wfs
end

function create_marker_instrs!(instr_lib, pulses, marker_chan)
	for p in pulses
		num_points = round(UInt64, length(p) * DAC_CLOCK)
		instr_lib[p] = Marker(marker_chan, num_points, p.amp > 0.5, true)
	end
end

function find_next_analog_entry!(entry, chan, wf_lib, analog_timestamps, id_ch, Z_chs_id)
	if any([id_ch[ct] > length(entry.pulses[chan[ct]]) for ct in 1:length(chan)])
		return
	end
	sim_chs_id = find(x-> x == minimum(analog_timestamps), analog_timestamps)
	filter!(x->!(x in Z_chs_id), sim_chs_id)
	pulses = []
	for ct in sim_chs_id
		pulse = entry.pulses[chan[ct]][id_ch[ct]]
		if typeof(pulse) == QGL.ZPulse
			if length(sim_chs_id) > 1 #if simultaneous logical channels, remove the Z pulses one by one
				push!(Z_chs_id, ct)
			end
			id_ch[ct]+=1
			return pulse, Z_chs_id
		end
		push!(pulses, pulse)
	end
	if length(sim_chs_id) == 1
		chan_select = sim_chs_id[1]
	else
		#select pulses on simultaneous channels
		nonid_ids = find([!wf_lib[pulse].isTA for pulse in pulses]) #find non-Id pulses
		if length(nonid_ids) > 1
			error("Only a single non-Id channel allowed")
		elseif length(nonid_ids) == 1
			chan_select = sim_chs_id[nonid_ids[1]]
		else #select the channel with the shortest Id
			chan_select = sim_chs_id[indmin([pulse.length for pulse in pulses])]
		end
	end
	next_entry = entry.pulses[chan[chan_select]][id_ch[chan_select]]
	for ct in sim_chs_id
		analog_timestamps[ct] += wf_lib[entry.pulses[chan[ct]][id_ch[ct]]].count + 1
		id_ch[ct]+=1
	end
	Z_chs_id = []
	return next_entry, Z_chs_id
end

function create_instrs(seqs, wf_lib, chans, chan_freqs)
	instrs = APS2Instruction[]

	# TODO: sort out whether we need any modulation commands
	# freqs = any(e.frequency != 0 for e in seqs if typeof(e) == QGL.Pulse)
	# frame_changes = any(typeof(e) == QGL.ZPulse for e in seqs)

	num_chans = length(chans)
	time_stamp = zeros(Int, num_chans)
	idx = ones(Int, num_chans)
	all_done = zeros(Bool, num_chans)
	num_entries = zeros(Int, num_chans)

	reset_phase_instr = modulation_instr(RESET_PHASE, 0x7)
	sync_instr = convert(APS2Instruction, QGL.sync())
	chan_freq_instrs = APS2Instruction[]
	if !isempty(chan_freqs)
		nco_select = Dict{QGL.Channel, UInt8}(chan => ct for (ct, chan) in enumerate(keys(chan_freqs)))
		for (chan, freq) in chan_freqs
			push!(chan_freq_instrs, modulation_instr(SET_FREQ, nco_select[chan], round(Int32, -freq / FPGA_CLOCK * 2^28 )) )
		end
	end

	for entry in seqs
		if typeof(entry) == QGL.PulseBlock
			# zero-out status vectors
			fill!(time_stamp, 0)
			fill!(idx, 1)
			for (ct, chan) in enumerate(chans)
				num_entries[ct] = minimum([length(entry.pulses[ch]) for ch in chan])
				all_done[ct] = num_entries[ct] == 0
			end

			analog_timestamps = zeros(Int, length(chan_freqs))
			analog_idx = ones(Int, length(chan_freqs))
			Z_chs_id = []
			# serialize pulses from the PulseBlock
			# round-robin through the channels until all are exhausted
			while !all(all_done)
				next_instr_time = minimum(time_stamp)

				for (ct, chan) in enumerate(chans)
					if (!all_done[ct]) && (time_stamp[ct] <= next_instr_time)

						if length(chan)>1 # multiple logical channels per analog channel
							next_entry_info = find_next_analog_entry!(entry, chan, wf_lib, analog_timestamps, analog_idx, Z_chs_id)
							if typeof(next_entry_info) == Void
								all_done[ct] = true
								break
							end
							next_entry, Z_chs_id = next_entry_info
						else
							next_entry = entry.pulses[chan[1]][idx[ct]]
							idx[ct] += 1
							all_done[ct] = idx[ct] > num_entries[ct]
						end
						if typeof(next_entry) == QGL.Pulse
							wf = wf_lib[next_entry]
							if typeof(next_entry.channel) == QGL.Qubit || typeof(next_entry.channel) == QGL.Edge
								# TODO: inject frequency update if necessary
								push!(instrs, modulation_instr(MODULATE, nco_select[next_entry.channel], wf.count))
							end
							push!(instrs, wf.instruction)
							time_stamp[ct] += wf.count+1
						elseif typeof(next_entry) == QGL.ZPulse
							# round phase to 28 bit integer
							fixed_pt_phase = round(Int32, mod(-next_entry.angle, 1) * 2^28 )
							push!(instrs, modulation_instr(UPDATE_FRAME, nco_select[next_entry.channel], fixed_pt_phase))
							# if the Z pulse is the first thing after a WAIT then we need to inject it before the WAIT
							if ((instrs[end-1] >> 60) & 0xff) == WAIT
								instrs[[end-1,end]] = instrs[[end, end-1]]
							end
						else
							error("Untranslated pulse block entry")
						end

					end

				end
			end

		else # if it's not a PulseBlock it must be a ControlFlow operation
			# convert control flow to APS2Instruction
			if entry.op == QGL.WAIT
				# heuristic to inject SYNC before a wait
				push!(instrs, sync_instr)
				# heuristic to reset modulation engine phase set frequency before wait for trigger
				push!(instrs, reset_phase_instr)
				for instr in chan_freq_instrs
					push!(instrs, instr)
				end
			end
			push!(instrs, convert(APS2Instruction, entry))
		end

	end

	return instrs

end

function convert(::Type{APS2Instruction}, cf::QGL.ControlFlow)
	if cf.op == QGL.WAIT
		return UInt64(WAIT << 4 | 0x1) << 56 | UInt64(WAIT_TRIG) << WFM_OP_OFFSET
	elseif cf.op == QGL.GOTO
		return UInt64(GOTO << 4 | 0x1) << 56 | UInt64(cf.target)
	elseif cf.op == QGL.SYNC
		return UInt64(SYNC << 4 | 0x1) << 56 | UInt64(WAIT_SYNC) << WFM_OP_OFFSET
	else
		error("Untranslated control flow instruction")
	end
end

function write_to_file(filename, instrs, wfs)
	# flatten waveforms to vector
	wf_vec = Vector{Complex{Int16}}()
	if !isempty(wfs)
		resize!(wf_vec, sum(length(wf) for wf in wfs))
		idx = 1
		for wf in wfs
			wf_vec[idx:idx+length(wf)-1] = wf
			idx += length(wf)
		end
	end

	#prepend the sequence directory
	seq_name_dir = filename[1:match(r"-", filename).offset-1]
	seq_path = joinpath(config.sequence_files_path, seq_name_dir)
	# create the sequence directory if necessary
	if !isdir(seq_path)
		mkpath(seq_path)
	end
	filename = joinpath(seq_path, filename)

	h5open(filename, "w") do f
		attrs(f)["Version"] = 4.0
		attrs(f)["target hardware"] = "APS2"
		attrs(f)["minimum firmware version"] = 4.0
		attrs(f)["channelDataFor"] = UInt16[1; 2]
		chan_1 = g_create(f, "chan_1")
		write(chan_1, "waveforms", real(wf_vec))
		write(chan_1, "instructions", instrs)
		chan_2 = g_create(f, "chan_2")
		write(chan_2, "waveforms", imag(wf_vec))
	end
end

end

ext: jl


text:
#---------------------------------------------------------
# # [Plots.jl examples](@id 2_Plots)
#---------------------------------------------------------

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/examples/2_Plots.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/examples/2_Plots.ipynb)

#md # !!! note
#md #     These examples are available as Jupyter notebooks.
#md #     You can execute them online with [binder](https://mybinder.org/) or just view them with [nbviewer](https://nbviewer.jupyter.org/) by clicking on the badges above!

# These examples were slightly modified from some of [the GR examples in the Plots.jl documentation](http://docs.juliaplots.org/latest/examples/gr/) and can be used as both a tutorial or as a series of test for the UnitfulRecipes package
# (they are essentially the same except we have added some units to the data).

# First we need to tell Julia we are using Plots, Unitful, and UnitfulRecipes

using Plots, Unitful, UnitfulRecipes

# ## Lines

plot(Plots.fakedata(50, 5) * u"m", w=3)

# ## Parametric plots

plot(t -> sin(t)*u"s", t -> sin(2t)*u"m", 0, 2π, line=4, leg=false, fill=(0, :orange))

# ## Colors

y = rand(100)*u"km"
plot((0:10:100)*u"hr", rand(11, 4)*u"km", lab="lines", w=3, palette=:grays, fill=0, α=0.6)
scatter!(y, zcolor=abs.(y .- 0.5u"km"), m=(:heat, 0.8, Plots.stroke(1, :green)), ms=10 * abs.(y .- 0.5u"km") .+ 4u"km", lab="grad")

# ## Global

# Note that a few changes had to be made for this to work.

using Statistics
y = rand(20, 3)*u"W"
x = (1:size(y,1))*u"Hz"
plot(x, y, xlabel="XLABEL", xlims=(-5, 30), xflip=true, xticks=0:2:20, background_color=RGB(0.2, 0.2, 0.2), leg=false)
hline!(mean(y, dims=1) + rand(1, 3)*u"W", line=(4, :dash, 0.6, [:lightgreen :green :darkgreen]))
vline!([5, 10]*u"Hz")
title!("TITLE")
yaxis!("YLABEL", :log10)

# ## Arguments

ys = Vector[rand(10), rand(20)] .* u"km"
plot(ys, color=[:black :orange], line=(:dot, 4), marker=([:hex :d], 12, 0.8, Plots.stroke(3, :gray)))

# ## Build plot in pieces

plot(rand(100) / 3 * u"km", reg=true, fill=(0, :green))
scatter!(rand(100) * u"km", markersize=6, c=:orange)

# ## Histogram2D

histogram2d(randn(10000) * u"cm", randn(10000) * u"cm", nbins=20)

# ## Line types

# ```
# linetypes = [:path :steppre :steppost :sticks :scatter]
# n = length(linetypes)
# x = Vector[sort(rand(20)) for i = 1:n] * u"km"
# y = rand(20, n) * u"ms"
# plot(x, y, line=(linetypes, 3), lab=map(string, linetypes), ms=15)
# ```

# ## Line styles

styles = intersect([:solid, :dash, :dot, :dashdot, :dashdotdot], Plots.supported_styles())
styles = reshape(styles, 1, length(styles))
n = length(styles)
y = cumsum(randn(20, n), dims=1) * u"km"
plot(y, line=(5, styles), label=map(string, styles), legendtitle="linestyle")

# ## Marker types

markers = intersect(Plots._shape_keys, Plots.supported_markers())
markers = reshape(markers, 1, length(markers))
n = length(markers)
x = (range(0, stop=10, length=n + 2))[2:end - 1] * u"km"
y = repeat(reshape(reverse(x), 1, :), n, 1)
scatter(x, y, m=(8, :auto), lab=map(string, markers), bg=:linen, xlim=(0, 10), ylim=(0, 10))

# ## Bar

bar(randn(99) * u"km")

# ## Histogram

histogram(randn(1000) * u"km", bins=:scott, weights=repeat(1:5, outer=200))

# ## Subplots

l = @layout([a{0.1h};b [c;d e]])
plot(randn(100, 5) * u"km", layout=l, t=[:line :histogram :scatter :steppre :bar], leg=false, ticks=nothing, border=:none)

# ## Adding to subplots

plot(Plots.fakedata(100, 10) * u"km", layout=4, palette=[:grays :blues :heat :lightrainbow], bg_inside=[:orange :pink :darkblue :black])

# ## Contour plots

x = (1:0.05:10) * u"m"
y = (1:0.01:2) * u"s"
f(x,y) = x^2 / y
z = f.(x',y)
p1 = contour(x, y, f, fill=true)
p2 = contour(x, y, z)
p3 = contourf(x, y, z)
plot(p1, p2, p3)

# ## 3D

n = 100
ts = range(0, stop=8π, length=n) * u"rad"
x = @. ts * cos(ts)
y = @. 0.1ts * sin(ts)
z = ts
plot(x, y, z, zcolor=reverse(z), m=(10, 0.8, :blues, Plots.stroke(0)), leg=false, cbar=true, w=5, xlabel="x", ylabel="y", zlabel="z")
plot!(zeros(n), zeros(n), z, w=5)

# ## Groups and Subplots

group = rand(map((i->"group $(i)"), 1:4), 100)
plot(rand(100)*u"km", layout=@layout([a b;c]), group=group, linetype=[:bar :scatter :steppre], linecolor=:match)

# ## Heatmap, categorical axes, and aspect_ratio

xs = [string("x", i) for i = 1:10]
ys = [string("y", i) for i = 1:4]
z = float((1:4) * reshape(1:10, 1, :)) * u"km"
heatmap(xs, ys, z, aspect_ratio=1)

# ## Magic grid argument

x = rand(10) * u"km"
p1 = plot(x, title="Default looks")
p2 = plot(x, grid=(:y, :olivedrab, :dot, 1, 0.9), title="Modified y grid")
p3 = plot(deepcopy(p2), title="Add x grid")
xgrid!(p3, :on, :cadetblue, 2, :dashdot, 0.4)
plot(p1, p2, p3, layout=(1, 3), label="", fillrange=0, fillalpha=0.3)

# ## Framestyle

# Suggestion: we might want to not add the unit label when the axis is not shown?

scatter(fill(randn(10), 6) * u"m", fill(randn(10), 6) * u"s", framestyle=[:box :semi :origin :zerolines :grid :none], title=[":box" ":semi" ":origin" ":zerolines" ":grid" ":none"], color=permutedims(1:6), layout=6, label="", markerstrokewidth=0, ticks=-2:2)

# ## Lines and markers with varying colors

# note that marker_z as a function did not work so it is modified here

t = range(0, stop=1, length=100) * u"s"
θ = 6π * u"rad/s" * t
x = @. t * cos(θ)
y = @. t * sin(θ)
z = x + y
p1 = plot(x, y, line_z=t, linewidth=3, legend=false)
p2 = scatter(x, y, marker_z=z, color=:bluesreds, legend=false)
plot(p1, p2)


text:
import Base:writemime

function html_body(p::VisJSGraph)
    """

    """
end

function get_color(c)
  if c != ""
    return ",color: \'$c\'"
  else
    return ""
  end
end

function generate_plot(p::VisJSGraph)
    nodes = ["{id:$(p.n[i].id),
               label: '$(p.n[i].label)',
               group: '$(p.n[i].group)'
               $(get_color(p.n[i].color))
              }" for i in 1:size(p.n,1)]
    edges = ""
    for item in p.e
      edges *= "{from:$(item[1]),to:$(item[2])},"
    end
    """require.config({
          paths: {
            vis: ["https://cdnjs.cloudflare.com/ajax/libs/vis/4.15.1/vis.min"],
          }
        });
    require(["vis"], function(vis){
    console.log(vis);
    // create an array with nodes
    var nodes = new vis.DataSet([
    $(join(nodes,","))
    ]);

    // create an array with edges
    var edges = new vis.DataSet([
    $(edges)
    ]);

    // create a network
    var container = document.getElementById('$(p.divid)');

    // provide the data in the vis format
    var data = {
        nodes: nodes,
        edges: edges
    };
    var options = {
      nodes: {
          shape: 'dot',
        },
    };

    // initialize your network!
    var network = new vis.Network(container, data, options);
});
    """
end

function Base.writemime(io::IO, ::MIME"text/html", p::VisJSGraph)
    print(io, html_body(p))
end


text:
module DocStringExtensionsMock

#If we don' thave DocStringExtensions installed, make a fall-back that prevents an error:
using Pkg
#If DocStringExtensions is installed, use it:
if Base.find_package("DocStringEstensions")!==nothing
    using DocStringExtensions
else
  TYPEDSIGNATURES=""
  SIGNATURES=""
end

export SIGNATURES, TYPEDSIGNATURES

end


text:
# This file is a part of BitOperations.jl, licensed under the MIT License (MIT).

# TODO: Deprecate custom byte-order array functions in favour of broadcasting.

function bswap! end
export bswap!

function ntoh! end
export ntoh!

function hton! end
export hton!

function ltoh! end
export ltoh!

function htol! end
export htol!

bswap!(x::AbstractArray) = begin
    @inbounds for i in eachindex(x)
        x[i] = Base.bswap(x[i])
    end
    x
end

bswap!(dest::AbstractArray, src) = begin
    @inbounds for i in eachindex(dest, src)
        dest[i] = Base.bswap(src[i])
    end
    dest
end

function _copyto_if_not_same!(dest::AbstractArray, src::AbstractArray)
    if dest !== src
        copyto!(dest, src)
    end
    dest
end

if ENDIAN_BOM == 0x01020304
    ntoh!(dest::AbstractArray, src) = _copyto_if_not_same!(dest, src)
    hton!(dest::AbstractArray, src) = _copyto_if_not_same!(dest, src)
    ltoh!(dest::AbstractArray, src) = bswap!(dest, src)
    htol!(dest::AbstractArray, src) = bswap!(dest, src)
elseif ENDIAN_BOM == 0x04030201
    ntoh!(dest::AbstractArray, src) = bswap!(dest, src)
    hton!(dest::AbstractArray, src) = bswap!(dest, src)
    ltoh!(dest::AbstractArray, src) = _copyto_if_not_same!(dest, src)
    htol!(dest::AbstractArray, src) = _copyto_if_not_same!(dest, src)
else
    error("Unknown machine endianess")
end


text:
module TestingMultivariateOnlineStatistics

using MultivariateOnlineStatistics, Test, Statistics, StatsBase

using MultivariateOnlineStatistics: storage

# Dummy private function used to "label" the tests.
check(flag::Bool, comment) = flag

@testset "MultivariateOnlineStatistics" begin
    T = Float64
    dims = (2, 3, 4)
    N = length(dims)
    numb = 17
    X = [rand(T, dims) for t in 1:numb]
    # Constructors and initialization
    A = IndependentStatistics{2,T}(dims)
    B = IndependentStatistics{2,T}(dims...)
    C = IndependentStatistics{2,T,length(dims)}(map(Int16, dims)...)
    avg = IndependentStatistics((Array{T}(undef, dims),))
    @test_throws ArgumentError IndependentStatistics{2,T,length(dims)+1}(dims)
    @test order(A) === 2
    @test order(avg) === 1
    @test eltype(A) === T
    @test ndims(A) === length(dims)
    @test size(A) === dims
    @test ntuple(k -> size(A, k), ndims(A)) === size(A)
    @test axes(A) === map(Base.OneTo, size(A))
    @test ntuple(k -> axes(A, k), ndims(A)) === axes(A)
    @test nobs(A) == 0
    @test extrema(mean(A)) === (zero(T), zero(T))
    @test extrema(std(A; corrected = false)) === (zero(T), zero(T))
    @test extrema(var(A; corrected = false)) === (zero(T), zero(T))
    @test_throws ErrorException std(A; corrected = true)
    @test_throws ErrorException var(A; corrected = true)
    push!(A, X[1])
    @test nobs(A) == 1
    @test mean(A) == X[1]
    @test extrema(std(A; corrected = false)) === (zero(T), zero(T))
    @test extrema(var(A; corrected = false)) === (zero(T), zero(T))
    @test_throws ErrorException std(A; corrected = true)
    @test_throws ErrorException var(A; corrected = true)
    push!(A, X[2])
    m1 = (X[1] + X[2])./2
    m2 = (X[1] - m1).^2 + (X[2] - m1).^2
    @test nobs(A) == 2
    @test mean(A) ≈ m1
    @test std(A; corrected = false) ≈ sqrt.(m2./2)
    @test var(A; corrected = false) ≈ m2./2
    @test std(A; corrected = true)  ≈ sqrt.(m2)
    @test var(A; corrected = true)  ≈ m2
    push!(A, X[3])
    m1 = (X[1] + X[2] + X[3])./3
    m2 = (X[1] - m1).^2 + (X[2] - m1).^2 + (X[3] - m1).^2
    @test nobs(A) == 3
    @test mean(A) ≈ m1
    @test std(A; corrected = false) ≈ sqrt.(m2./3)
    @test var(A; corrected = false) ≈ m2./3
    @test std(A; corrected = true)  ≈ sqrt.(m2./2)
    @test var(A; corrected = true)  ≈ m2./2
    push!(A, X[4])
    m1 = (X[1] + X[2] + X[3] + X[4])./4
    m2 = (X[1] - m1).^2 + (X[2] - m1).^2 + (X[3] - m1).^2 + (X[4] - m1).^2
    @test nobs(A) == 4
    @test mean(A) ≈ m1
    @test std(A; corrected = false) ≈ sqrt.(m2./4)
    @test var(A; corrected = false) ≈ m2./4
    @test std(A; corrected = true)  ≈ sqrt.(m2./3)
    @test var(A; corrected = true)  ≈ m2./3
    # Check element-by-element mean.
    arr = mean(A)
    flag = true
    for i in eachindex(arr)
        flag &= (arr[i] ≈ mean(A, i))
    end
    @test check(flag, "mean(A, i::Integer) == mean(A)[i]")
    flag = true
    for i in CartesianIndices(arr)
        flag &= (arr[i] ≈ mean(A, i))
    end
    @test check(flag, "mean(A, i::CartesianIndex) == var(A)[i]")
    flag = true
    for i in CartesianIndices(arr)
        flag &= (arr[i] ≈ mean(A, Tuple(i)))
    end
    @test check(flag, "mean(A, i::Tuple{Vararg{Integer}}) == var(A)[i]")
    # Check element-by-element variance.
    arr = var(A)
    flag = true
    for i in eachindex(arr)
        flag &= (arr[i] ≈ var(A, i))
    end
    @test check(flag, "var(A, i::Integer) == var(A)[i]")
    flag = true
    for i in CartesianIndices(arr)
        flag &= (arr[i] ≈ var(A, i))
    end
    @test check(flag, "var(A, i::CartesianIndex) == var(A)[i]")
    flag = true
    for i in CartesianIndices(arr)
        flag &= (arr[i] ≈ var(A, Tuple(i)))
    end
    @test check(flag, "var(A, i::Tuple{Vararg{Integer}}) == var(A)[i]")
    arr = var(A; corrected=false)
    flag = true
    for i in eachindex(arr)
        flag &= (arr[i] ≈ var(A, i; corrected=false))
    end
    @test check(flag, "var(A, i::Integer; corrected=false) == var(A; corrected=false)[i]")
    # Check element-by-element standard deviation.
    arr = std(A)
    flag = true
    for i in eachindex(arr)
        flag &= (arr[i] ≈ std(A, i))
    end
    @test check(flag, "std(A, i::Integer) == std(A)[i]")
    flag = true
    for i in CartesianIndices(arr)
        flag &= (arr[i] ≈ std(A, i))
    end
    @test check(flag, "std(A, i::CartesianIndex) == std(A)[i]")
    flag = true
    for i in CartesianIndices(arr)
        flag &= (arr[i] ≈ std(A, Tuple(i)))
    end
    @test check(flag, "std(A, i::Tuple{Vararg{Integer}}) == std(A)[i]")
    arr = std(A; corrected=false)
    flag = true
    for i in eachindex(arr)
        flag &= (arr[i] ≈ std(A, i; corrected=false))
    end
    @test check(flag, "std(A, i::Integer; corrected=false) == std(A; corrected=false)[i]")

    # Push remaining terms in two different ways.
    push!(A, X[5:7]...)
    merge!(A, X[8:end])
    @test nobs(A) == length(X)
    # Push all terms at once.
    merge!(B, X)
    @test nobs(B) == length(X)
    for (a, b) in zip(storage(A), storage(B))
        @test a == b
    end
    merge!(avg, X)
    @test nobs(avg) == length(X)
    @test mean(avg) == mean(B)
    # Push in reverse order.
    for i in reverse(1:length(X))
       push!(C, X[i])
    end
    @test nobs(C) == nobs(A)
    for (a, c) in zip(storage(A), storage(C))
        @test a ≈ c
    end
    # Empty statistics.
    empty!(B)
    @test eltype(B) === T
    @test ndims(B) === length(dims)
    @test nobs(B) == 0
    for b in storage(B)
        @test extrema(b) === (zero(T), zero(T))
    end
    empty!(C)
    # Collect statistics by parts.
    n1 = length(X)÷3
    n2 = length(X)÷2
    for i in 1:n1
        push!(B, X[i])
    end
    for i in n1+1:n1+n2
        push!(C, X[i])
    end
    for i in n1+n2+1:length(X)
        push!(B, X[i])
    end
    merge!(B, C)
    @test nobs(B) == nobs(A)
    for (a, b) in zip(storage(A), storage(B))
        @test a ≈ b
    end
    # Check merge! into an empty statistics.
    empty!(B)
    merge!(B, A)
    @test nobs(B) == nobs(A)
    for (a, b) in zip(storage(A), storage(B))
        @test a == b
    end
    # Check merge! from a single sample statistics.
    empty!(A)
    empty!(B)
    for x in X
        push!(A, x)
        merge!(B, push!(empty!(C), x))
        if nobs(A) ≥ 5
            break
        end
    end
    @test nobs(B) == nobs(A)
    for (a, b) in zip(storage(A), storage(B))
        @test a == b
    end

    # Partially qualified constructors.
    merge!(empty!(A), X)
    for B in (IndependentStatistics(storage(A), nobs(A)),
              IndependentStatistics{order(A)}(storage(A), nobs(A)),
              IndependentStatistics{order(A),eltype(A)}(storage(A), nobs(A)),
              IndependentStatistics{order(A),eltype(A),ndims(A)}(storage(A), nobs(A)),)
        @test nobs(B) == length(X)
        @test order(B) == order(A)
        @test eltype(B) == eltype(A)
        @test size(B) == size(A)
    end
    for B in (IndependentStatistics(storage(A)),
              IndependentStatistics{order(A)}(storage(A)),
              IndependentStatistics{order(A),eltype(A)}(storage(A)),
              IndependentStatistics{order(A),eltype(A),ndims(A)}(storage(A)),)
        @test nobs(B) == 0
        @test order(B) == order(A)
        @test eltype(B) == eltype(A)
        @test size(B) == size(A)
    end
end

end # module

max_issues_repo_name: emmt/MultivariateOnlineStatistics

max_forks_repo_name: emmt/MultivariateOnlineStatistics


text:
module SpringCollab2020GlassBerryModule

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/glassBerry" GlassBerry(x::Integer, y::Integer, checkpointID::Integer=-1, order::Integer=-1, nodes::Array{Tuple{Integer, Integer}, 1}=Tuple{Integer, Integer}[])

const placements = Ahorn.PlacementDict(
	"Glass Strawberry (Spring Collab 2020)" => Ahorn.EntityPlacement(
		GlassBerry,
		"point"
	),
)

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

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

	nodes = get(entity.data, "nodes", ())
	hasPips = length(nodes) > 0

	sprite = "collectables/SpringCollab2020/glassBerry/idle00"
	seedSprite = "collectables/strawberry/seed00"

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

	for node in nodes
		nx, ny = node

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

	return res
end

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

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

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

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

	nodes = get(entity.data, "nodes", ())
	hasPips = length(nodes) > 0

	sprite = "collectables/SpringCollab2020/glassBerry/idle00"
	seedSprite = "collectables/strawberry/seed00"

	for node in nodes
		nx, ny = node

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

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

end


text:
using Documenter, JuMPTutorials

makedocs(
    modules = [JuMPTutorials],
    format = Documenter.HTML(),
    sitename = "JuMPTutorials.jl",
    authors  = "Arpit Bhatia",
    pages = [
        "Home" => "index.md",
        "Adding a New Tutorial" => "new.md",
        "Notes" => "notes.md",
        "Future Work" => "future.md",
        "Function Index" => "api.md"
    ]
)

Documenter.deploydocs(repo = "github.com/JuliaOpt/JuMPTutorials.jl.git")


text:
################################################################################
#
#   DEFINITION OF 3D UNITCELLS
#
#   FOLLOWING UNITCELLS ARE IMPLEMENTED
#   - cubic
#
#   Code is included in subfiles
#
################################################################################

# Cubic lattice
include("unitcells_3d/cubic.jl")

# FCC lattice (face centered cubic)
include("unitcells_3d/fcc.jl")
# diamond lattice
include("unitcells_3d/diamond.jl")
# Pyrochlore lattice
include("unitcells_3d/pyrochlore.jl")

# Pyrochlore lattice
include("unitcells_3d/hyperkagome.jl")

# (8,3)a lattice
include("unitcells_3d/8_3_a.jl")
# (8,3)b lattice
include("unitcells_3d/8_3_b.jl")
# (8,3)c lattice
include("unitcells_3d/8_3_c.jl")
# (8,3)n lattice
include("unitcells_3d/8_3_n.jl")

# (9,3)a lattice
include("unitcells_3d/9_3_a.jl")

# (10,3)a lattice
include("unitcells_3d/10_3_a.jl")
# (10,3)b lattice
include("unitcells_3d/10_3_b.jl")
# (10,3)c lattice
include("unitcells_3d/10_3_c.jl")
# (10,3)d lattice
include("unitcells_3d/10_3_d.jl")


text:
using KrylovMethods
using Base.Test
using LinearOperators
include("getDivGrad.jl")

println("=== Testing CG ===")
# small full system
A = [4.0 1; 1 4]
rhs = [2.0;2]
x,flag,relres,iter,resvec = cg(A,rhs,tol=1e-15,out=2,storeInterm=true)
@test norm(A*x[:,end]-rhs)/norm(rhs) <= 1e-15

# test message and flag when stopping early
x,flag,relres,iter,resvec = cg(A,randn(2),tol=1e-15,out=2,maxIter=1)
@test flag==-1

# test for zero right hand-side
x,flag = cg(A,zeros(2),out=2)
@test flag == -9
@test all(x.==0)

# negative definite matrix
b = [1 -1.0]
KKT = [A b'; b 0]
x,flag, = cg(KKT,randn(3),out=1)
@test flag==-2

# CG: test sparse Laplacian
A = getDivGrad(8,8,8)
Alinop = LinearOperator(A)
L = tril(A)
D = diag(A)
U = triu(A)
n = size(A,1)
JAC = (x) -> D.\x
SGS = (x) -> L\(D.*(U\x))

rhs = randn(size(A,1))
tolCG = 1e-5
# tests with A being matrix
xCG,flagCG,relresCG,iterCG,resvecCG       = cg(A,rhs,tol=tolCG,maxIter=100)
xJAC,flagJAC,relresJAC,iterJAC,resvecJAC  = cg(A,rhs,tol=tolCG,maxIter=100,M=JAC,out=1)
xSGS,flagSGS,relresSGS,iterSGS,resvecSGS  = cg(A,rhs,tol=tolCG,maxIter=100,M=SGS)
# tests with A being function
xCGmf,flagCG,relresCG,iterCG,resvecCG       = cg(Alinop,rhs,tol=tolCG,maxIter=100)
xJACmf,flagJAC,relresJAC,iterJAC,resvecJAC  = cg(Alinop,rhs,tol=tolCG,maxIter=100,M=JAC)
xSGSmf,flagSGS,relresSGS,iterSGS,resvecSGS  = cg(Alinop,rhs,tol=tolCG,maxIter=100,M=SGS)
# tests with random starting guess
xCGr,flagCGr,relresCGr,iterCGr,resvecCGr       = cg(Alinop,rhs,tol=tolCG,maxIter=100,x=randn(size(rhs)))
xJACr,flagJACr,relresJACr,iterJACr,resvecJACr  = cg(Alinop,rhs,tol=tolCG,maxIter=100,M=JAC,x=randn(size(rhs)))
xSGSr,flagSGSr,relresSGSr,iterSGSr,resvecSGSr  = cg(A,rhs,tol=tolCG,maxIter=100,x=randn(size(rhs)))

# test relative residuals
@test norm(A*xCG-rhs)/norm(rhs) <= tolCG
@test norm(A*xSGS-rhs)/norm(rhs) <= tolCG
@test norm(A*xJAC-rhs)/norm(rhs) <= tolCG
@test norm(A*xCGmf-rhs)/norm(rhs) <= tolCG
@test norm(A*xSGSmf-rhs)/norm(rhs) <= tolCG
@test norm(A*xJACmf-rhs)/norm(rhs) <= tolCG
@test norm(A*xCGr-rhs)/norm(rhs) <= tolCG*10
@test norm(A*xSGSr-rhs)/norm(rhs) <= tolCG*10
@test norm(A*xJACr-rhs)/norm(rhs) <= tolCG*10
# preconditioners should at least not increase number of iter
@test iterJAC==iterCG
@test iterSGS<=iterJAC

println("=== CG : All tests passed. ====")


text:
# From StatsMakie
function _linear(x::AbstractVector{T}, y::AbstractVector;
                 npoints = 100, wts = similar(x, 0),
                 interval = length(wts) > 0 ? nothing : :confidence) where T
    # Note: confidence intervals are currently not supported for WLS in GLM.jl
    try
        y = collect(y)
        x = collect(x)
        wts = collect(wts)
        lin_model = GLM.lm([ones(T, length(x)) x], y, wts=wts)
        x_min, x_max = extrema(x)
        x_new = range(x_min, x_max, length = npoints)
        pred = GLM.predict(lin_model,
                           [ones(T, npoints) x_new],
                           interval=interval)
        if !isnothing(interval)
            y_new, lower, upper = pred
        else
            y_new = pred
        end
        # the GLM predictions always return matrices
        x, y = x_new, vec(y_new)
        lines = mapping(x, y)
        if !isnothing(interval)
            band = mapping(x, vec(lower), vec(upper))
            return visual(:Lines) * lines + visual(:Band, alpha = 0.2) * band
        else
            return visual(:Lines) * lines
        end
    catch e
        @warn "Linear fit not possible for the given data"
        return AlgebraicList()
    end
end

"""
    linear(x, y; wts = similar(x, 0), interval = length(wts) > 0 ? nothing : :confidence)

Compute a linear fit of `y ~ 1 + x`. Weighted data is supported via the keyword `wts`.
Use `interval` to specify what type of interval the shaded band should represent.
Valid values of interval are `:confidence` delimiting the uncertainty of the predicted
relationship, and `:prediction` delimiting estimated bounds for new data points.
"""
const linear = Analysis(_linear)

# TODO: multidimensional case as a heatmap or surface plot
function _smooth(x, y; npoints = 100, kwargs...)
    min, max = extrema(x)
    min < max || return AlgebraicList()
    model = Loess.loess(Float64.(x), Float64.(y); kwargs...)
    us = collect(range(min, stop = max, length = npoints))
    vs = Loess.predict(model, us)
    return visual(:Lines) * mapping(us, vs)
end

"""
    smooth(x, y, span=0.75, degreee=2)

Fit a loess model. `span` is the degree of smoothing, typically in `[0,1]`.
Smaller values result in smaller local context in fitting. `degree` is the polynomial
degree used in the loess model.
"""
const smooth = Analysis(_smooth)


text:
module Train

export train!

using FluxGAN: GAN
using FluxGAN: image_grid

using Flux
using Flux.Losses: logitbinarycrossentropy
using EllipsisNotation
using Statistics
using Colors
using FileIO
using CUDA

######################
# training functions #
######################

# main training function - attempts to train on GPU

function train!(model::GAN, train_tensor::AbstractArray; kws...)
    if has_cuda_gpu()
        println("training GAN with GPU...")
        train_gpu!(model, train_tensor; kws...)
    else
        println("no GPU available :( training GAN with CPU...")
        train_cpu!(model, train_tensor; kws...)
    end
end

# gpu training function

function train_gpu!(model::GAN, train_tensor::AbstractArray;
                    iterations=1000, verbose=true, skip=50,
                    gif=false, gif_path="gifs", gif_layout=(5, 5), gif_fps=25,
                    gif_filename="", return_losses=false, kws...)

    d = model.hparams.latent_dim
    m = model.hparams.minibatch

    mean_gen_losses = Float32[]
    mean_dscr_losses = Float32[]

    gen_losses = Float32[]
    dscr_losses = Float32[]

    t0 = time_ns()

    if gif
        rows, cols = gif_layout
        h, w = model.hparams.img_size[1:2]
        gif_tensor = Array{Color}(undef, h * rows, w * cols, div(iterations, skip))
    end

    model = fmap(gpu, model)
    for i = 1:iterations
        if verbose && i % skip == 0
            mean_gen_loss = mean(gen_losses)
            mean_dscr_loss = mean(dscr_losses)
            push!(mean_gen_losses, mean_gen_loss)
            push!(mean_dscr_losses, mean_dscr_loss)
            training_message(i, iterations, t0, mean_gen_loss, mean_dscr_loss, skip; kws...)
            if gif
                gif_tensor[:,:,div(i, skip)] = image_grid(model, gif_layout; seed=true)
            end
            gen_losses = []
            dscr_losses = []
            t0 = time_ns()
        end
        for _ = 1:model.hparams.dscr_loops
            x_tensor = train_tensor[.., rand(1:size(train_tensor)[end], m)] |> gpu
            z_tensor = CUDA.randn(Float32, d, m)
            dscr_loss = train_discriminator!(model, x_tensor, z_tensor)
            push!(dscr_losses, dscr_loss)
        end
        z_tensor = CUDA.randn(Float32, d, m)
        gen_loss = train_generator!(model, z_tensor)
        push!(gen_losses, gen_loss)
    end
    model = fmap(cpu, model)

    if gif
        save(gif_path * "/" * gif_filename * "_n_$(iterations)_" *
             "grid_$(rows)_$(cols)" * ".gif",
             gif_tensor; fps=gif_fps)
    end

    if return_losses
        return mean_gen_losses, mean_dscr_losses
    end
end

# cpu training function

function train_cpu!(model::GAN, train_tensor::AbstractArray;
                    iterations=1000, verbose=true, skip=50,
                    gif=false, gif_path="images/gifs", gif_layout=(5, 5), gif_fps=25,
                    gif_filename="", kws...)

    d = model.hparams.latent_dim
    m = model.hparams.minibatch

    gen_loss = 0
    dscr_loss = 0

    t0 = time_ns()

    if gif
        gif_tensor = Array{Color}(undef, size(train_tensor)[1:2]...,
                                       *(gif_layout...),
                                       div(iterations, skip))
    end

    for i = 1:iterations
        if verbose && i % skip == 0
            training_message(i, iterations, t0, gen_loss, dscr_loss, skip; kws...)
            if gif
                gif_tensor[..,div(i, skip)] = image_grid_tensor(model,
                                                *(gif_layout...);
                                                seed=true)
            end
            gen_loss = 0
            dscr_loss = 0
            t0 = time_ns()
        end
        for _ = 1:model.hparams.dscr_loops
            x_tensor = train_tensor[.., rand(1:size(train_tensor)[end], m)]
            z_tensor = randn(Float32, d, m)
            dscr_loss += train_discriminator!(model, x_tensor, z_tensor)
        end
        z_tensor = randn(Float32, d, m)
        gen_loss += train_generator!(model, z_tensor)
    end

    if gif
        save(gif_path * "/" * gif_filename * "_n_$(iterations)_" *
             "grid_$(gif_layout[1])_$(gif_layout[2])" * ".gif",
             gif_tensor; fps=gif_fps)
    end
end

# discriminator train function

function train_discriminator!(model::GAN, xs::AbstractArray, zs::AbstractArray)
    ω = Flux.params(model.D)
    loss, ∇ = Flux.pullback(ω) do
        discriminator_loss(model.D(xs), model.D(model.G(zs)))
    end
    ∇ = ∇(1f0)
    Flux.update!(model.D_opt, ω, ∇)
    return loss
end

# generator train function

function train_generator!(model::GAN, zs::AbstractArray)
    θ = Flux.params(model.G)
    loss, ∇ = Flux.pullback(θ) do
        generator_loss(model.D(model.G(zs)))
    end
    ∇ = ∇(1f0)
    Flux.update!(model.G_opt, θ, ∇)
    return loss
end

##################
# loss functions #
##################

function discriminator_loss(real_outputs, fake_outputs)
    real_loss = logitbinarycrossentropy(real_outputs, 1)
    fake_loss = logitbinarycrossentropy(fake_outputs, 0)
    return real_loss + fake_loss
end

generator_loss(fake_outputs) = logitbinarycrossentropy(fake_outputs, 1)

# training message function

# TODO: add plotting functionality
#           * mean training loss
#           * maybe error bars

function training_message(i, iterations, t0, gen_loss, dscr_loss, skip)
    t = time_ns()
    println("\niteration: $i of ", iterations)
    println("avg seconds per loop   = ", Float32((t - t0) / skip)*1f-9)
    println("avg generator loss     = ", Float32(gen_loss))
    println("avg discriminator loss = ", Float32(dscr_loss))
end

end

max_stars_repo_name: aarontrowbridge/GAN

max_issues_repo_name: aarontrowbridge/GAN

max_forks_repo_name: aarontrowbridge/GAN


text:
module ChromeProfileFormat

import JSON, Profile

include("core.jl")

end # module


text:
using Test
using MPI

MPI.Init()

comm_size = MPI.Comm_size(MPI.COMM_WORLD)
comm_rank = MPI.Comm_rank(MPI.COMM_WORLD)

# send to next higher process, with wraparound
dest = mod(comm_rank+1, comm_size)
src  = mod(comm_rank-1, comm_size)

# test simple type
mutable struct NotABits
    a::Any
end

@testset "Non bitstype" begin
    @test_throws ArgumentError MPI.Datatype(NotABits)
end

struct Boundary
    c::UInt16  # force some padding to be inserted
    a::Int
    b::UInt8
end
@testset "Compound type" begin
    sz = sizeof(Boundary)
    al = Base.datatype_alignment(Boundary)
    @test MPI.Types.extent(MPI.Datatype(Boundary)) == (0, cld(sz,al)*al)

    arr = [Boundary( (comm_rank + i) % 127, i + comm_rank, i % 64) for i = 1:3]
    req_send = MPI.Isend(arr, dest, 1, MPI.COMM_WORLD)

    # receive the message
    arr_recv = Array{Boundary}(undef, 3)
    req_recv = MPI.Irecv!(arr_recv, src, 1, MPI.COMM_WORLD)

    MPI.Wait!(req_send)
    MPI.Wait!(req_recv)

    # check received array
    for i=1:3
        bndry_i = arr_recv[i]
        @test bndry_i.a == (src + i)
        @test bndry_i.b == i % 64
        @test bndry_i.c == (src + i) % 127
    end
end

struct Boundary2
    a::UInt32
    b::Tuple{Int, UInt8}
    c::Nothing
end
@testset "nested types" begin
    sz = sizeof(Boundary2)
    al = Base.datatype_alignment(Boundary2)
    @test MPI.Types.extent(MPI.Datatype(Boundary2)) == (0, cld(sz,al)*al)

    arr = [Boundary2( (comm_rank + i) % 127, ( Int(i + comm_rank), UInt8(i % 64)), nothing) for i = 1:3]
    arr_recv = Array{Boundary2}(undef,3)

    req_send = MPI.Isend(arr, dest, 1, MPI.COMM_WORLD)
    req_recv = MPI.Irecv!(arr_recv, src, 1, MPI.COMM_WORLD)

    MPI.Wait!(req_send)
    MPI.Wait!(req_recv)

    # check received array
    for i=1:3
        bndry_i = arr_recv[i]
        @test bndry_i.a == (src + i) % 127
        @test bndry_i.b[1] == (src + i)
        @test bndry_i.b[2] == (i % 64)
        @test bndry_i.c === nothing
    end
end

primitive type Primitive16 16 end
primitive type Primitive24 24 end
primitive type Primitive80 80 end

@testset for PrimitiveType in (Primitive16, Primitive24, Primitive80)
    sz = sizeof(PrimitiveType)
    al = Base.datatype_alignment(PrimitiveType)
    @test MPI.Types.extent(MPI.Datatype(PrimitiveType)) == (0, cld(sz,al)*al)

    if VERSION < v"1.3" && PrimitiveType == Primitive80
        # alignment is broken on earlier Julia versions
        continue
    end

    arr = [Core.Intrinsics.trunc_int(PrimitiveType, UInt128(comm_rank + i)) for i = 1:4]
    arr_recv = Array{PrimitiveType}(undef,4)

    recv_req = MPI.Irecv!(arr_recv, src, 2, MPI.COMM_WORLD)
    send_req = MPI.Isend(arr, dest, 2, MPI.COMM_WORLD)

    MPI.Wait!(recv_req)
    MPI.Wait!(send_req)

    @test arr_recv == [Core.Intrinsics.trunc_int(PrimitiveType, UInt128(src + i)) for i = 1:4]
end

@testset "packed non-aligned tuples" begin
    T = NTuple{3,UInt8}

    sz = sizeof(T)
    al = Base.datatype_alignment(T)
    @test MPI.Types.extent(MPI.Datatype(T)) == (0, cld(sz,al)*al)

    arr = [(UInt8(comm_rank),UInt8(i),UInt8(0)) for i = 1:8]
    arr_recv = Array{T}(undef,8)

    req_send = MPI.Isend(arr, dest, 1, MPI.COMM_WORLD)
    req_recv = MPI.Irecv!(arr_recv, src, 1, MPI.COMM_WORLD)

    MPI.Wait!(req_send)
    MPI.Wait!(req_recv)

    # check received array
    @test arr_recv == [(UInt8(src),UInt8(i),UInt8(0)) for i = 1:8]
end

@testset "0-sized type" begin
    sz = sizeof(Nothing)
    al = Base.datatype_alignment(Nothing)

    # OpenMPI gives incorrect values
    # see https://github.com/open-mpi/ompi/issues/7266
    # @test MPI.Types.extent(MPI.Datatype(Nothing)) == (0, cld(sz,al)*al)

    arr = [nothing for i = 1:100]
    arr_recv = Array{Nothing}(undef,100)

    req_send = MPI.Isend(arr, dest, 1, MPI.COMM_WORLD)
    req_recv = MPI.Irecv!(arr_recv, src, 1, MPI.COMM_WORLD)

    MPI.Wait!(req_send)
    MPI.Wait!(req_recv)

    # check received array
    @test arr_recv == [nothing for i = 1:100]
end

MPI.Barrier(MPI.COMM_WORLD)
MPI.Finalize()
@test MPI.Finalized()


text:
function [X,J,dXdx]=JacobiAnnulus(ksi,eta,F,Grid)
n1=size(ksi,1);
n2=size(eta,2);
X=zeros(n1,n2,2);
Phi1=atan2(F.P(2,1),F.P(1,1));
Rad1=sqrt(F.P(1,1)^2+F.P(2,1)^2);
Phi2=atan2(F.P(2,2),F.P(1,2));
Rad2=sqrt(F.P(1,2)^2+F.P(2,2)^2);
Phi3=atan2(F.P(2,3),F.P(1,3));
Rad3=sqrt(F.P(1,3)^2+F.P(2,3)^2);
Phi4=atan2(F.P(2,4),F.P(1,4));
Rad4=sqrt(F.P(1,4)^2+F.P(2,4)^2);
if Phi1<0
  Phi1=Phi1+2*pi;
end
if Phi2<0
  Phi2=Phi2+2*pi;
end
if Phi3<0
  Phi3=Phi3+2*pi;
end
if Phi4<0
  Phi4=Phi4+2*pi;
end
if Phi1>Phi3
  Phi3=Phi3+2*pi;
end
if Phi2>Phi4
  Phi4=Phi4+2*pi;
end
Phi=0.25*(Phi1*(1-ksi)*(1-eta)+...
  Phi2*(1+ksi)*(1-eta)+...
  Phi3*(1+ksi)*(1+eta)+...
  Phi4*(1-ksi)*(1+eta));
Rad=0.25*(Rad1*(1-ksi)*(1-eta)+...
  Rad2*(1+ksi)*(1-eta)+...
  Rad3*(1+ksi)*(1+eta)+...
  Rad4*(1-ksi)*(1+eta));
X(:,:,1)=cos(Phi).*Rad ;
X(:,:,2)=sin(Phi).*Rad;

if nargout > 2
  dX1dPhi=-sin(Phi).*Rad;
  dX1dRad=cos(Phi);
  dX2dPhi=cos(Phi).*Rad;
  dX2dRad=sin(Phi);

  dPhidksi=0.25*((-Phi1+Phi2+...
    Phi3-Phi4)*ones(n1,1)*ones(1,n2)+...
    ( Phi1-Phi2+...
    Phi3-Phi4)*ones(n1,1)*eta);
  dPhideta=0.25*((-Phi1-Phi2+...
    Phi3+Phi4)*ones(n1,1)*ones(1,n2)+...
    ( Phi1-Phi2+...
    Phi3-Phi4)*ksi*ones(1,n2));
  dRaddksi=0.25*((-Rad1+Rad2+...
    Rad3-Rad4)*ones(n1,1)*ones(1,n2)+...
    ( Rad1-Rad2+...
    Rad3-Rad4)*ones(n1,1)*eta);
  dRaddeta=0.25*((-Rad1-Rad2+...
    Rad3+Rad4)*ones(n1,1)*ones(1,n2)+...
    ( Rad1-Rad2+...
    Rad3-Rad4)*ksi*ones(1,n2));

  dXdx=zeros(n1,n2,2,2);
  dXdx(:,:,1,1)=dX1dPhi.*dPhidksi+dX1dRad.*dRaddksi;
  dXdx(:,:,2,1)=dX2dPhi.*dPhidksi+dX2dRad.*dRaddksi;
  dXdx(:,:,1,2)=dX1dPhi.*dPhideta+dX1dRad.*dRaddeta;
  dXdx(:,:,2,2)=dX2dPhi.*dPhideta+dX2dRad.*dRaddeta;
  J=dXdx(:,:,1,1).*dXdx(:,:,2,2)-dXdx(:,:,1,2).*dXdx(:,:,2,1);
end
end


text:
module QSpin

using SparseArrays

import Base: kron, +, -, *

include("spinstate.jl")

export SpinState
export tensorprod

end

max_stars_repo_name: amiragha/QSpin

max_issues_repo_name: amiragha/QSpin

max_forks_repo_name: amiragha/QSpin


text:
"""
    TwoDVort(x,z,s,vc,vx,vz)

Defines a 2D vortex at `(x,z)` with vorticity `s` and vortex-core
radius `vc`.

`vx` and `vz` are induced velocity at centroid used in vortex
interaction calculations

"""
mutable struct  TwoDVort
    x :: Float64
    z :: Float64
    s :: Float64
    vc :: Float64
    vx :: Float64
    vz :: Float64
end

struct TwoDFlowField
    velX :: MotionDef
    velZ :: MotionDef
    u :: Vector{Float64}
    w :: Vector{Float64}
    tev :: Vector{TwoDVort}
    lev :: Vector{TwoDVort}
    extv :: Vector{TwoDVort}
    function TwoDFlowField(velX = ConstDef(0.), velZ = ConstDef(0.))
        u = [0;]
        w = [0;]
        tev = TwoDVort[]
        lev = TwoDVort[]
        extv = TwoDVort[]
        new(velX, velZ, u, w, tev, lev, extv)
    end
end

struct TwoDSurf
    c :: Float64
    uref :: Float64
    coord_file :: String
    pvt :: Float64
    ndiv :: Int16
    naterm :: Int16
    kindef :: KinemDef
    cam :: Vector{Float64}
    cam_slope :: Vector{Float64}
    theta :: Vector{Float64}
    x :: Vector{Float64}
    kinem :: KinemPar
    bnd_x :: Vector{Float64}
    bnd_z :: Vector{Float64}
    uind :: Vector{Float64}
    wind :: Vector{Float64}
    downwash :: Vector{Float64}
    a0 :: Vector{Float64}
    aterm :: Vector{Float64}
    a0dot :: Vector{Float64}
    adot :: Vector{Float64}
    a0prev :: Vector{Float64}
    aprev :: Vector{Float64}
    bv :: Vector{TwoDVort}
    lespcrit :: Vector{Float64}
    levflag :: Vector{Int8}
    initpos :: Vector{Float64}
    rho :: Float64

    function TwoDSurf(coord_file, pvt, kindef,lespcrit=zeros(1); c=1., uref=1., ndiv=70, naterm=35, initpos = [0.; 0.], rho = 0.04)
        theta = zeros(ndiv)
        x = zeros(ndiv)
        cam = zeros(ndiv)
        cam_slope = zeros(ndiv)
        bnd_x = zeros(ndiv)
        bnd_z = zeros(ndiv)
        kinem = KinemPar(0, 0, 0, 0, 0, 0)

        dtheta = pi/(ndiv-1)
        for ib = 1:ndiv
            theta[ib] = (ib-1.)*dtheta
            x[ib] = c/2. *(1-cos(theta[ib]))
        end
        if (coord_file != "FlatPlate")
            cam, cam_slope = camber_calc(x, coord_file)
        end
        #camspl = Spline1D(theta, cam)
        #for i = 1:ndiv
        #    cam_slope[i] = Dierckx.derivative(camspl, theta[i])
        #end

        if (typeof(kindef.alpha) == EldUpDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c
        elseif (typeof(kindef.alpha) == EldRampReturnDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c
        elseif (typeof(kindef.alpha) == ConstDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = 0.
        elseif (typeof(kindef.alpha) == SinDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c
        elseif (typeof(kindef.alpha) == CosDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c
        end

        if (typeof(kindef.h) == EldUpDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        elseif (typeof(kindef.h) == EldUpIntDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        elseif (typeof(kindef.h) == EldRampReturnDef)
            kinem.h= kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        elseif (typeof(kindef.h) == ConstDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = 0.
        elseif (typeof(kindef.h) == SinDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        elseif (typeof(kindef.h) == CosDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        end

        if (typeof(kindef.u) == EldUpDef)
            kinem.u = kindef.u(0.)*uref
            kinem.udot = ForwardDiff.derivative(kindef.u,0.)*uref*uref/c
        elseif (typeof(kindef.u) == EldRampReturnDef)
            kinem.u, kinem.udot = kindef.u(0.)
            kinem.u = kinem.u*uref
            kinem.udot = kinem.udot*uref*uref/c
        elseif (typeof(kindef.u) == ConstDef)
            kinem.u = kindef.u(0.)*uref
            kinem.udot = 0.
        end

        for i = 1:ndiv
            bnd_x[i] = -((c - pvt*c)+((pvt*c - x[i])*cos(kinem.alpha))) + (cam[i]*sin(kinem.alpha)) + initpos[1]
            bnd_z[i] = kinem.h + ((pvt*c - x[i])*sin(kinem.alpha))+(cam[i]*cos(kinem.alpha)) + initpos[2]
        end
        uind = zeros(ndiv)
        wind = zeros(ndiv)
        downwash = zeros(ndiv)
        a0 = zeros(1)
        a0dot = zeros(1)
        aterm = zeros(naterm)
        adot = zeros(naterm)
        a0prev = zeros(1)
        aprev = zeros(naterm)
        bv = TwoDVort[]
        for i = 1:ndiv-1
            push!(bv,TwoDVort(0,0,0,0.02*c,0,0))
        end
        levflag = [0]
        new(c, uref, coord_file, pvt, ndiv, naterm, kindef, cam, cam_slope, theta, x, kinem, bnd_x, bnd_z, uind, wind, downwash, a0, aterm, a0dot, adot, a0prev, aprev, bv, lespcrit, levflag, initpos, rho)
    end
end

mutable struct TwoDOFPar
    x_alpha :: Float64
    r_alpha :: Float64
    kappa :: Float64
    w_alpha :: Float64
    w_h :: Float64
    w_alphadot :: Float64
    w_hdot :: Float64
    cubic_h_1 :: Float64
    cubic_h_3 :: Float64
    cubic_alpha_1 :: Float64
    cubic_alpha_3 :: Float64
end

#For airfoil with 2DOF in pitch and plunge
mutable struct KinemPar2DOF
    alpha :: Float64
    h :: Float64
    alphadot :: Float64
    hdot :: Float64
    u :: Float64
    udot :: Float64
    alphaddot :: Float64
    hddot :: Float64
    alpha_pr :: Float64
    h_pr :: Float64
    alphadot_pr :: Float64
    alphadot_pr2 :: Float64
    alphadot_pr3 :: Float64
    hdot_pr :: Float64
    hdot_pr2 :: Float64
    hdot_pr3 :: Float64
    alphaddot_pr :: Float64
    alphaddot_pr2 :: Float64
    alphaddot_pr3 :: Float64
    hddot_pr :: Float64
    hddot_pr2 :: Float64
    hddot_pr3 :: Float64

    function KinemPar2DOF(alpha, h, alphadot, hdot, u)
        new(alpha, h, alphadot, hdot, u, 0., 0., 0., alpha, h, alphadot, alphadot, alphadot, hdot, hdot, hdot, 0., 0., 0., 0., 0., 0.)
    end
end

struct KelvinCondition
    surf :: TwoDSurf
    field :: TwoDFlowField
end

function (kelv::KelvinCondition)(tev_iter::Array{Float64})

    nlev = length(kelv.field.lev)
    ntev = length(kelv.field.tev)

    uprev, wprev = ind_vel([kelv.field.tev[ntev]], kelv.surf.bnd_x, kelv.surf.bnd_z)

    #Update the TEV strength
    kelv.field.tev[ntev].s = tev_iter[1]

    unow, wnow = ind_vel([kelv.field.tev[ntev]], kelv.surf.bnd_x, kelv.surf.bnd_z)

    kelv.surf.uind[:] = kelv.surf.uind[:] .- uprev .+ unow
    kelv.surf.wind[:] = kelv.surf.wind[:] .- wprev .+ wnow

    #Calculate downwash
    update_downwash(kelv.surf, [kelv.field.u[1],kelv.field.w[1]])

    #Calculate first two fourier coefficients
    update_a0anda1(kelv.surf)

    val = kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0[1] + kelv.surf.aterm[1]/2.) -
        kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0prev[1] + kelv.surf.aprev[1]/2.) +
        kelv.field.tev[ntev].s

    return val
end

struct KelvinConditionMult
    surf :: Vector{TwoDSurf}
    field :: TwoDFlowField
end

function (kelv::KelvinConditionMult)(tev_iter::Array{Float64})

    nsurf = length(kelv.surf)

    val = zeros(nsurf)

    nlev = length(kelv.field.lev)
    ntev = length(kelv.field.tev)

    tev_list = kelv.field.tev[ntev-nsurf+1:ntev]
    for i = 1:nsurf

        bv_list = TwoDVort[]
        for j = 1:nsurf
            if i != j
                bv_list = [bv_list; kelv.surf[j].bv]
            end
        end
        uprev, wprev = ind_vel([tev_list; bv_list], kelv.surf[i].bnd_x, kelv.surf[i].bnd_z)

        #Update the TEV strength
        kelv.field.tev[ntev-nsurf+i].s = tev_iter[i]

        unow, wnow = ind_vel([tev_list; bv_list], kelv.surf[i].bnd_x, kelv.surf[i].bnd_z)

        kelv.surf[i].uind[:] = kelv.surf[i].uind[:] .- uprev .+ unow
        kelv.surf[i].wind[:] = kelv.surf[i].wind[:] .- wprev .+ wnow

        #Calculate downwash
        update_downwash(kelv.surf[i], [kelv.field.u[1],kelv.field.w[1]])
    end

    for i = 1:nsurf
        #Update Fourier coefficients and bv strength
        update_a0anda1(kelv.surf[i])
        update_a2toan(kelv.surf[i])
        update_bv(kelv.surf[i])

        val[i] = kelv.surf[i].uref*kelv.surf[i].c*pi*(kelv.surf[i].a0[1] + kelv.surf[i].aterm[1]/2.) -
            kelv.surf[i].uref*kelv.surf[i].c*pi*(kelv.surf[i].a0prev[1] + kelv.surf[i].aprev[1]/2.) +
            kelv.field.tev[ntev-nsurf+i].s
    end
    return val
end

struct KelvinKutta
    surf :: TwoDSurf
    field :: TwoDFlowField
end

function (kelv::KelvinKutta)(v_iter::Array{Float64})
    val = zeros(2)

    nlev = length(kelv.field.lev)
    ntev = length(kelv.field.tev)

    uprev, wprev = ind_vel([kelv.field.tev[ntev]; kelv.field.lev[nlev]], kelv.surf.bnd_x, kelv.surf.bnd_z)

    #Update the TEV and LEV strengths
    kelv.field.tev[ntev].s = v_iter[1]
    kelv.field.lev[nlev].s = v_iter[2]

    unow, wnow = ind_vel([kelv.field.tev[ntev]; kelv.field.lev[nlev]], kelv.surf.bnd_x, kelv.surf.bnd_z)

    kelv.surf.uind[:] = kelv.surf.uind[:] .- uprev .+ unow
    kelv.surf.wind[:] = kelv.surf.wind[:] .- wprev .+ wnow

    #Calculate downwash
    update_downwash(kelv.surf ,[kelv.field.u[1],kelv.field.w[1]])

    #Calculate first two fourier coefficients
    update_a0anda1(kelv.surf)

    val[1] = kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0[1] + kelv.surf.aterm[1]/2.) -
        kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0prev[1] + kelv.surf.aprev[1]/2.) +
        kelv.field.tev[ntev].s + kelv.field.lev[nlev].s

    if (kelv.surf.a0[1] > 0)
        lesp_cond = kelv.surf.lespcrit[1]
    else
        lesp_cond = -kelv.surf.lespcrit[1]
    end
    val[2] = kelv.surf.a0[1]-lesp_cond

    return val
end

struct KelvinKuttaMult
    surf :: Vector{TwoDSurf}
    field :: TwoDFlowField
    shed_ind :: Vector{Int}
end

function (kelv::KelvinKuttaMult)(v_iter::Array{Float64})
    nsurf = length(kelv.surf)
    nshed = length(kelv.shed_ind)
    val = zeros(nsurf+nshed)

    nlev = length(kelv.field.lev)
    ntev = length(kelv.field.tev)

    tev_list = kelv.field.tev[ntev-nsurf+1:ntev]
    lev_list = kelv.field.lev[nlev-nsurf.+kelv.shed_ind]

    levcount = 0
    for i = 1:nsurf
        bv_list = TwoDVort[]
        for j = 1:nsurf
            if i != j
                bv_list = [bv_list; kelv.surf[j].bv]
            end
        end
        uprev, wprev = ind_vel([tev_list; lev_list; bv_list], kelv.surf[i].bnd_x, kelv.surf[i].bnd_z)

        #Update the shed vortex strengths
        kelv.field.tev[ntev-nsurf+i].s = v_iter[i]

        if i in kelv.shed_ind
            levcount += 1
            kelv.field.lev[nlev-nsurf+i].s = v_iter[nsurf+levcount]
        end

        unow, wnow = ind_vel([tev_list; lev_list; bv_list], kelv.surf[i].bnd_x, kelv.surf[i].bnd_z)

        kelv.surf[i].uind[:] = kelv.surf[i].uind[:] .- uprev .+ unow
        kelv.surf[i].wind[:] = kelv.surf[i].wind[:] .- wprev .+ wnow

        #Calculate downwash
        update_downwash(kelv.surf[i], [kelv.field.u[1],kelv.field.w[1]])
    end

    levcount = 0
    for i = 1:nsurf
        #Update Fourier coefficients and bv strength
        update_a0anda1(kelv.surf[i])
        update_a2toan(kelv.surf[i])
        update_bv(kelv.surf[i])

        if i in kelv.shed_ind
            val[i] = kelv.surf[i].uref*kelv.surf[i].c*pi*(kelv.surf[i].a0[1] + kelv.surf[i].aterm[1]/2.) -
                kelv.surf[i].uref*kelv.surf[i].c*pi*(kelv.surf[i].a0prev[1] + kelv.surf[i].aprev[1]/2.) +
                kelv.field.tev[ntev-nsurf+i].s + kelv.field.lev[nlev-nsurf+i].s
            levcount += 1
            if (kelv.surf[i].a0[1] > 0)
                lesp_cond = kelv.surf[i].lespcrit[1]
            else
                lesp_cond = -kelv.surf[i].lespcrit[1]
            end

            val[levcount+nsurf] = kelv.surf[i].a0[1] - lesp_cond
        else
            val[i] = kelv.surf[i].uref*kelv.surf[i].c*pi*(kelv.surf[i].a0[1] + kelv.surf[i].aterm[1]/2.) -
                kelv.surf[i].uref*kelv.surf[i].c*pi*(kelv.surf[i].a0prev[1] + kelv.surf[i].aprev[1]/2.) +
                kelv.field.tev[ntev-nsurf+i].s
        end
    end
    return val
end

max_stars_repo_name: hjabird/UNSflow

max_issues_repo_name: hjabird/UNSflow

max_forks_repo_name: hjabird/UNSflow


text:
"""
    RSOCBridge{T, F, G}

The `RotatedSecondOrderCone` is `SecondOrderCone` representable; see [1, p. 104].
Indeed, we have ``2tu = (t/√2 + u/√2)^2 - (t/√2 - u/√2)^2`` hence
```math
2tu \\ge || x ||_2^2
```
is equivalent to
```math
(t/√2 + u/√2)^2 \\ge || x ||_2^2 + (t/√2 - u/√2)^2.
```
We can therefore use the transformation ``(t, u, x) \\mapsto (t/√2+u/√2, t/√2-u/√2, x)``.
Note that the linear transformation is a symmetric involution (i.e. it is its own transpose and its own inverse).
That means in particular that the norm is of constraint primal and duals are preserved by the tranformation.

[1] Ben-Tal, Aharon, and Arkadi Nemirovski. *Lectures on modern convex optimization: analysis, algorithms, and engineering applications*. Society for Industrial and Applied Mathematics, 2001.
"""
struct RSOCBridge{T, F, G} <: SetMapBridge{T, MOI.SecondOrderCone, MOI.RotatedSecondOrderCone, F, G}
    constraint::CI{F, MOI.SecondOrderCone}
end

function rotate_function_type(G::Type{<:MOI.AbstractVectorFunction}, T::Type)
    S = MOIU.promote_operation(/, T, MOIU.scalar_type(G), T)
    Y = MOIU.promote_operation(-, T, S, S)
    Z = MOIU.promote_operation(+, T, S, S)
    return MOIU.promote_operation(vcat, T, Z, Y, G)
end

function concrete_bridge_type(::Type{<:RSOCBridge{T}},
                              G::Type{<:MOI.AbstractVectorFunction},
                              ::Type{MOI.RotatedSecondOrderCone}) where T
    return RSOCBridge{T, rotate_function_type(G, T), G}
end

function map_set(::Type{<:RSOCBridge}, set::MOI.RotatedSecondOrderCone)
    return MOI.SecondOrderCone(MOI.dimension(set))
end
function inverse_map_set(::Type{<:RSOCBridge}, set::MOI.SecondOrderCone)
    return MOI.RotatedSecondOrderCone(MOI.dimension(set))
end

"""
    SOCRBridge{T, F, G}

We simply do the inverse transformation of [`RSOCBridge`](@ref). In fact, as the
transformation is an involution, we do the same transformation.
"""
struct SOCRBridge{T, F, G} <: SetMapBridge{T, MOI.RotatedSecondOrderCone, MOI.SecondOrderCone, F, G}
    constraint::CI{F, MOI.RotatedSecondOrderCone}
end

function concrete_bridge_type(::Type{<:SOCRBridge{T}},
                              G::Type{<:MOI.AbstractVectorFunction},
                              ::Type{MOI.SecondOrderCone}) where T
    return SOCRBridge{T, rotate_function_type(G, T), G}
end

function map_set(::Type{<:SOCRBridge}, set::MOI.SecondOrderCone)
    return MOI.RotatedSecondOrderCone(MOI.dimension(set))
end
function inverse_map_set(::Type{<:SOCRBridge}, set::MOI.RotatedSecondOrderCone)
    return MOI.SecondOrderCone(MOI.dimension(set))
end

function map_function(::Type{<:Union{RSOCBridge{T}, SOCRBridge{T}}}, func) where T
    scalars = MOIU.eachscalar(func)
    t = scalars[1]
    u = scalars[2]
    x = scalars[3:end]
    s2 = √T(2)
    ts = MOIU.operate!(/, T, t, s2)
    us = MOIU.operate!(/, T, u, s2)
    # Cannot use `operate!` here since `ts` and `us` are needed for the next
    # line
    y  = ts - us
    z  = MOIU.operate!(+, T, ts, us)
    return MOIU.operate(vcat, T, z, y, x)
end
# The map is an involution
inverse_map_function(BT::Type{<:Union{RSOCBridge, SOCRBridge}}, func) = map_function(BT, func)
# The map is symmetric
adjoint_map_function(BT::Type{<:Union{RSOCBridge, SOCRBridge}}, func) = map_function(BT, func)
# The map is a symmetric involution
inverse_adjoint_map_function(BT::Type{<:Union{RSOCBridge, SOCRBridge}}, func) = map_function(BT, func)


text:
#could be written more efficiently
#but is not meant to be evaluated repeatedly like the likelihood
#so not a priority

function viterbi(model::DynamicDiscreteModel,data::Array{Int,1})
	T=length(data)
	dx=size(model.mu,1)

	memory=Array(Int,dx,T-1)

	model.psi[:]=model.mu[:,data[1]]/sum(model.mu[:,data[1]])
	for t=2:T
		mm=reshape(model.m[:,data[t-1],:,data[t]],(dx,dx))
		for jx=1:dx
			model.psi[jx],memory[jx,t-1]=findmax(model.psi.*mm[:,jx])
		end
		model.psi[:]=model.psi/sum(model.psi)
	end
	viterbipath=Array(Int,T)
	viterbipath[T]=indmax(model.psi)
	for t=T-1:-1:1
		viterbipath[t]=memory[viterbipath[t+1],t]
	end
	viterbipath
end

#wrapper for panel data
function viterbi(model::DynamicDiscreteModel,data::Array{Array,1})
	n=length(data)
	viterbipaths=Array(Array{Int,1},n)
	for i=1:n
		viterbipaths[i]=viterbi(model,data[i])
	end
	viterbipaths
end


text:
using MacroTools

Base.getindex(n::Nullable) = n.value

abstract AR

const metadata = MetaData()

immutable TableMap
    table::Table
    columns::Dict{Symbol, String}
    relations::Dict{Symbol, Dict{Symbol, Any}}
end

const table_map = Dict{Any, TableMap}()

function columnname(tbl, sqlname)
    for (k,v) in table_map[tbl].columns
        if v==sqlname
            return k
        end
    end
end

@enum State Transient Pending Persistent Detached
@enum ARSetState Unloaded Loaded

type ARSet{T<:AR}
    R::Vector{T}
    state::Symbol
    parent::AR
    column::Symbol
    state::ARSetState
    ARSet(parent, column) = new(T[], :unloaded, parent, column, Unloaded)
end

function Base.getindex(s::ARSet, idx)
    s.R[idx]
end

Base.endof(s::ARSet) = endof(s.R)
Base.length(s::ARSet) = length(s.R)

function Base.push!(s::ARSet, r::AR)
    push!(s.R, r)
    parent = s.parent
    if parent._state ∈ (Pending, Persistent)
        push!(parent._session[], r)
    end
    if parent._state == Persistent
        setforeign!(r, parent, s.column)
    end
    s
end

macro table(desc)
    success = @capture desc begin
        type T_
            tablename_
            fields__
        end
    end
    success || error("Invalid input to @table: $desc")
    new_fields = []
    columns = Column[]
    keys=[]
    name_map = Dict{Symbol, String}()
    field_types=Dict()
    relations = Dict{Symbol, Dict{Symbol, Any}}()
    for field in fields
        res = @match field begin
            (key_::(id_, t_, kwargs__)) => (:hasname, key, id, t, kwargs)
            (key_::(t_, kwargs__)) => (:noname, key, t, kwargs)
            (key_::relation(t_, kwargs__)) => (:relation, key, t, kwargs)
        end
        kind, res = res[1], res[2:end]
        if kind == :relation
            key, foreign_t, kwargs = res
            push!(new_fields, :($key::ARSet{$(esc(foreign_t))}))
            push!(keys, (:foreign, key, foreign_t))
            relations[key] = Dict{Symbol, Any}()
            for arg in kwargs
                @capture arg begin
                    a_ = b_
                end
                relations[key][a]=b
            end
        else
            if kind == :hasname
                key, id, t, kwargs = res
            elseif kind == :noname
                key, t, kwargs = res
                id = String(key)
            else
                error("Field $field not understood")
            end
            name_map[key] = id
            field_types[key] = convert(JuliaType, eval(t))
            # field_name = field_types[key].name.name
            field_name = field_types[key].name.name
            push!(new_fields, :($key::Nullable{$field_name}))

            push!(keys, (:field, key))
            kw = Dict()
            for arg in kwargs
                @capture arg begin
                    a_ = b_
                end
                kw[a] =b
            end
            push!(columns, Column(id, eval(t); kw...))
        end
    end
    defaults=[]
    for (idx, val) in enumerate(keys)
        kind = val[1]
        val = val[2:end]
        if kind == :field
            key = val[1]
            push!(defaults, :(this.$key=Nullable{$(field_types[key])}()))
        elseif kind == :foreign
            key, t = val
            push!(defaults, :(this.$key=ARSet{$(esc(t))}(this, $(QuoteNode(key)))))
        end
    end

    quote
        type $(esc(T)) <: AR
            $(new_fields...)
            _state::State
            _session::Nullable{Session}
            _key::Nullable{Int}

            function $(esc(T))(;kwargs...)
                this = new()
                $(defaults...)
                for (key, val) in kwargs
                    this.(key) = Nullable(val)
                end
                this._state = Transient
                this._session = Nullable{Session}()
                this._key = Nullable{Int}()
                this
            end
        end

        # function Base.show(io::IO, t::$(esc(T)))
        #     print(io, $(QuoteNode(T)), ": ")
        #     names = ["$key=>$(field_show(t[key]))" for key in $keys]
        #     print(io, join(names, ", "))
        # end

        SQLAlchemy.table_map[$(esc(T))] =
            SQLAlchemy.TableMap(SQLAlchemy.Table($tablename, metadata, ($columns)...),
                                $name_map, $relations)
    end
end

immutable ID
    key::Int
    table::DataType
end

function primarykey{T<:AR}(a::Type{T})
    tbl = table_map[a]
    pkeys = inspect(tbl.table)[:primary_key][:columns][:keys]()
    length(pkeys) == 0 && error("Wrapped class $(a) must have a primary key.")
    length(pkeys) > 1 && error("Composite primary keys not supported. Error on table $(a).")
    pkey = first(pkeys)
    for (k,v) in tbl.columns
        if v==pkey
            return k
        end
    end
    error("Internal error")
end

macro typelift(f)
    quote
        $(esc(f))(x::AR) = $f(typeof(x))
    end
end

@typelift primarykey

function ID{T<:AR}(t::T)
    pkeyname = primarykey(T)
    ID(t[pkeyname][], T)
end

type Session
    db::Connection
    new_::Set{AR}
    dirty::Set{AR}
    deleted::Set{AR}
    identity_map::Dict{ID, AR}
end

Base.in{T<:AR}(ar::T, s::Session) = ID(ar) ∈ s
Base.in(id::ID, s::Session) = id ∈ s.identity_map

function setstate!(a::AR, state::State)
    a._state = state
    if state == Persistent
        relations = table_map[typeof(a)].relations
        for relation in keys(relations)
            s = a[relation]
            for r in s.R
                setforeign!(r, a, s.column)
            end
        end
    end
    a
end

function dirty!(r::AR)
    if !isnull(r._session)
        push!(r._session[].dirty, r)
    end
end

function joincolumns(child_t, parent_t, column)
    table = gettable(child_t)
    parent_table = gettable(parent_t)
    columns = table.o[:columns][:items]()
    for (name, column) in columns
        for key in column[:foreign_keys]
            if key[:column][:table] == parent_table.o
                return columnname(child_t, name), columnname(parent_table, key[:column][:name])
            end
        end
    end
end

function setforeign!(child, parent, column)
    child_key, parent_key = joincolumns(typeof(child), typeof(parent), column)
    child[child_key] = parent[parent_key]
    child
end

function foreignkeys{T<:AR}(::Type{T})
    table_map[T].relations |> keys
end
@typelift foreignkeys

function Base.push!(s::Session, a::AR)
    a._state == Persistent && a ∈ s && return s
    a._state == Pending && a ∈ s.new_ && return s
    setstate!(a, Pending)
    push!(s.new_, a)
    a._session = s
    for keyname in foreignkeys(a)
        for r in a[keyname].R
            push!(s, r)
        end
    end
    s
end

function ensureloaded!{T}(session::Session, a::ARSet{T})
    a.state == Loaded && return a
    parent_t = typeof(a.parent)
    child_key, parent_key = joincolumns(T, parent_t, :a)
    a.R = query(session, T)|>filter(parent_t[parent_key] == T[child_key])|>all
    a.state = Loaded
    a
end

function Base.getindex(a::AR, field)
    val = a.(field)
    if a._state == Persistent && isa(val, ARSet)
        ensureloaded!(a._session[], val)
    end
    val
end

function Base.setindex!(a::AR, value, field)
    setindex!(a, Nullable(value), field)
end

function Base.setindex!(a::AR, value::Nullable, field)
    a.(field) = value
    dirty!(a)
end

Session(engine::Engine) = Session(connect(engine), Set{AR}(), Set{AR}(),
                                  Set{AR}(), Dict{ID, AR}())

function Base.show(io::IO, s::Session)
    print(io, "Session")
end

type SessionQuery
    session::Session
    table::Type
    select::Select
end

function query{T<:AR}(session::Session, ::Type{T})
    SessionQuery(session, T, select([gettable(T)]))
end

function Base.join{T<:AR}(::Type{T})
    query->begin
        t1 = table_map[query.table].table
        t2 = table_map[T].table
        clause = query.select |> selectfrom(join(t1, t2))
        SessionQuery(query.session, query.table, clause)
    end
end

function Base.all(s::SessionQuery)
    flush(s.session)
    recordset = s.session.db(s.select) |> fetchall
    ars = Vector{s.table}()
    for record in recordset.records
        push!(ars, push!(s.session, record, s.table))
    end
    ars
end

function Base.first(s::SessionQuery)
    flush(s.session)
    record = s.session.db(s.select) |> fetchone
    push!(s.session, record, s.table)
end

function fields(r::AR)
    column_map = table_map[typeof(r)].columns
    d = Dict()
    for (k, v) in column_map
        d[v] = r[k]
    end
    d
end

function sqlfields(r::AR)
    f = Dict()
    for (k, v) in fields(r)
        k = Symbol(k)
        if isnull(v)
            continue
        else
            f[k] = v[]
        end
    end
    f
end

gettable{T<:AR}(::Type{T}) = table_map[T].table
@typelift gettable

function Base.flush(s::Session)
    db = s.db
    new_ = [s.new_...]
    # for r in s.new_
    while !isempty(s.new_)
        r = pop!(s.new_)
        tbl = gettable(r)
        f = sqlfields(r)
        res = db(insert(tbl); f...)
        pkeys = res[:inserted_primary_key]
        isempty(pkeys) && error("Need primary key on $r")
        pkey = first(pkeys)
        r[primarykey(typeof(r))] = pkey
        s.identity_map[ID(r)] = r
        setstate!(r, Persistent)
    end
    empty!(s.new_)

    for r in s.dirty
        tbl = gettable(r)
        pkey = primarykey(r)
        stmt = update(tbl) |> where(tbl[pkey]==r[pkey][])
        db(stmt; sqlfields(r)...)
    end
    empty!(s.dirty)
end

Base.close(s::Session) = close(s.db)

function update!(a::AR, r::Record)
    column_map = table_map[typeof(a)].columns
    for (k,v) in column_map
        a[k] = r[Symbol(v)]
    end
    a
end

function Base.push!(s::Session, r::Record, tbl)
    pname = primarykey(tbl)
    sqlname = Symbol(table_map[tbl].columns[pname])
    id = ID(r[sqlname][], tbl)
    if id ∈ keys(s.identity_map)
        a = s.identity_map[id]
    else
        a = tbl()
        a._state = Persistent
        s.identity_map[id] = a
    end
    update!(a, r)
end

function Base.getindex{T<:SQLAlchemy.AR}(::Type{T}, jl_field)
    sql_field = table_map[T].columns[jl_field]
    gettable(T)[field]
end

function Base.filter(clause::Wrapped)
    query->begin
        new_clause = query.select |> where(clause)
        SessionQuery(query.session, query.table, new_clause)
    end
end

function orderby(q::SessionQuery, field)
    SessionQuery(q.session, q.table, q.select |> orderby(field))
end

# function Base.show{T<:AR}(io::IO, ::Type{T})
#     table = table_map[T].table
#     show(io, table)
# end


text:
using LinearAlgebra
using BlockFactorizations
using BenchmarkTools

d = 512
n, m = 16, 8
A = [randn(d, d) for i in 1:n, j in 1:m]
B = BlockFactorization(A)

x = [randn(d) for _ in 1:m]
y = [zeros(d) for _ in 1:n]
@btime A*x
@btime mul!(y, A, x);

x = randn(d*m)
y = zeros(d*n)
@btime B*x;
@btime mul!(y, B, x);

d = 512
n, m = 16, 8
# testing with diagonal components
A = [Diagonal(randn(d)) for i in 1:n, j in 1:m]
B = BlockFactorization(A)

x = [randn(d) for _ in 1:m]
y = [zeros(d) for _ in 1:n]
@btime A*x
@btime mul!(y, A, x);

x = randn(d*m)
y = zeros(d*n)
@btime B*x;
@btime mul!(y, B, x);


text:
function perform_dijkstra_fm(W, pstart; method="fm", svg_rate=10, niter=Inf, bound="sym")
    """
        # perform_fm_dijstra - *slow* (matlab) implementation of Dijstra and FM
        #
        #   [D,Dsvg,Ssvg] = perform_fm_dijstra(W, pstart, options);
        #
        #   W is an (n,n) metric matrix.
        #   pstart is a (2,k) starting points.
        #   options.method is either 'fm' or 'dijstra'
        #
        #   D is the final distance map to pstart
        #   options.svg_rate gives the rate at wich Dsvg and Ssvg is filled.
        #   options.niter can be used to limit the total number of steps (partial propagation).
        #
        #   Copyright (c) 2012 Gabriel Peyre
    """

##
# Size.

n = size(W,1)

##
# The four displacement vector to go to the four neightbors.

neigh = [[1;0] [-1;0] [0;1] [0;-1]]

##
# For simplicity of implementation, we use periodic boundary conditions.

function symmetrize(x,n)
    if (x<=0)
        x = 2-x;
    elseif (x>n)
        x = 2*n-x
    end
    return x
end

if bound == "per"
    boundary = x -> mod(x-1,n)+1
elseif bound == "sym"
    boundary = x -> [symmetrize(x[1],n), symmetrize(x[2],n)]
else
    error("Works only for per and sym.")
end
#     case 'sym'
#         boundary = @(x)x.*(x<=n & x>0) + (2-x).*(x<=0) + (2*n-x).*(x>n);
#     otherwise
#         error('Works only for per and sym.');
# end

##
# For a given grid index |k|, and a given neighboring index k in \({1,2,3,4}\),
# |Neigh(k,i)| gives the corresponding grid neigboring index.

ind2sub1 = k -> Int.([rem(k-1, n)+1; (k - rem(k-1, n) - 1)/n + 1])
sub2ind1 = u -> Int((u[2]-1)*n + u[1])
Neigh = (k,i) -> sub2ind1(boundary(ind2sub1(k) + neigh[:,i]));
extract   = (x,I) -> x[I]
extract1d = (x,I) -> extract(vec(x),I);

##
# Stack of starting points.
nstart = size(pstart,2)
I = Int.(zeros(nstart))
for i=1:nstart
    I[i] = sub2ind1(pstart[:,i])
end

##
# Initialize the distance to \(+\infty\), excepted for the boundary conditions.

D = zeros(n,n)+Inf; # current distance
for i in 1:nstart
    D[pstart[1,i], pstart[2,i]] = 0;
end

##
# Initialize the state to 0 (unexplored), excepted for the boundary point to \(1\)
# (front).

S = zeros(n,n);
for i in 1:nstart
    S[pstart[1,i],pstart[2,i]] = 1 # open
end

##
# Run!
compt=0
iter = 0
q = 100  # maximum number of saves
Dsvg = zeros(n,n,q)
Ssvg = zeros(n,n,q)
while !(isempty(I)) && iter<=niter
    iter = iter+1
    # pop from stack
    j = sortperm(extract1d(D,I))
    j=j[1]
    i = I[j]
    deleteat!(I,j)
    # declare dead
    u = ind2sub1(i);
    S[u[1],u[2]] = -1;
    # Make a list of neighbors that are not dead
    J = []
    for k in 1:4
        j = Neigh(i,k)
        if extract1d(S,j)!= -1
            # add to the list of point to update
            append!(J,j)
            if extract1d(S,j)==0
                # add to the front
                u = ind2sub1(j)
                S[u[1],u[2]] = 1
                append!(I,j)
            end
        end
    end
    # update neighbor values
    DNeigh = (D,k) -> extract1d(D,Neigh(j,k))
    for j in J
        dx = min(DNeigh(D,1), DNeigh(D,2))
        dy = min(DNeigh(D,3), DNeigh(D,4))
        u = ind2sub1(j)
        w = extract1d(W,j);
        if method=="dijkstra"
            D[u[1],u[2]] = min(dx + w, dy + w)
        else
            Delta = 2*w - (dx-dy).^2
            if (Delta>=0)
                D[u[1],u[2]] = (dx + dy + sqrt(Delta))/ 2
            else
                D[u[1],u[2]] = min(dx + w, dy + w)
            end
        end
    end
    # svd
    t = Base.div(iter,svg_rate)
    if mod(iter,svg_rate)==0 && (t<q)
        Dsvg[:,:,t] = D
        Ssvg[:,:,t] = S
    end
end
return (D,Dsvg,Ssvg)
end


text:
using Test

using HydrophoneCalibrations
using HydrophoneCalibrations: enforce_hertz_unit, interp
using HydrophoneCalibrations.Unitful

tolerance_v2p = 1e4u"Pa/V";

@testset "HydrophoneCalibrations.jl" begin
    @testset "precision acoustics hydrophone" begin
        v2p = volt_to_pressure(1e6, :PA_3197)
        sensitivity_at_1mhz = 34u"mV/MPa"
        conversion_factor = uconvert(u"Pa/V", 1/sensitivity_at_1mhz)
        @test isapprox(v2p, conversion_factor, atol=tolerance_v2p)
    end

    @testset "onda hydrophone" begin
        v2p = volt_to_pressure(1e6, :Onda_HGL0200_2322)
        sensitivity_at_1mhz = 4.11e-8u"V/Pa"
        conversion_factor = uconvert(u"Pa/V", 1/sensitivity_at_1mhz)
        @test isapprox(v2p, conversion_factor, atol=tolerance_v2p)
    end

    @testset "onda hydrophone with preamplifier (linear calculations)" begin
        v2p = volt_to_pressure(1e6, :Onda_HGL0200_2322, :Onda_AH2020_1238_20dB)
        sensitivity_at_1mhz = 4.11e-8u"V/Pa"
        hydrophone_capacitance_at_1mhz = 9.963u"pF"
        preamp_capacitance_at_1mhz = 7.08u"pF"
        preamp_gain_at_1mhz = 19.92u"dB"
        capacitance_factor = hydrophone_capacitance_at_1mhz/(hydrophone_capacitance_at_1mhz+preamp_capacitance_at_1mhz)
        system_sensitivity = uconvertrp(NoUnits, preamp_gain_at_1mhz) * sensitivity_at_1mhz * capacitance_factor
        conversion_factor = uconvert(u"Pa/V", 1/system_sensitivity)
        @test isapprox(v2p, conversion_factor, atol=tolerance_v2p)
    end

    @testset "onda hydrophone with preamplifier (logarithmic calculations)" begin
        v2p = volt_to_pressure(1e6, :Onda_HGL0200_2322, :Onda_AH2020_1238_20dB)
        sensitivity_at_1mhz = 4.11e-8u"V/Pa"
        hydrophone_capacitance_at_1mhz = 9.963u"pF"
        preamp_capacitance_at_1mhz = 7.08u"pF"
        preamp_gain_at_1mhz = 19.92u"dB"
        capacitance_factor = hydrophone_capacitance_at_1mhz/(hydrophone_capacitance_at_1mhz+preamp_capacitance_at_1mhz)
        system_sensitivity = preamp_gain_at_1mhz + ustrip(uconvertrp(u"dBV/Pa", sensitivity_at_1mhz))*u"dBV" + uconvertrp(u"dB", capacitance_factor)
        system_sensitivity = uconvertrp(u"V", system_sensitivity) / u"Pa"
        conversion_factor = uconvert(u"Pa/V", 1/system_sensitivity)
        @test isapprox(v2p, conversion_factor, atol=tolerance_v2p)
    end

    @testset "onda hydrophone with unknown preamplifier" begin
        @test_throws MethodError volt_to_pressure(1e6, :Onda_HGL0200_2322, :my_dummy)
    end

    @testset "volt_to_pressure_and_phase() with unknown preamp" begin
        @test_throws MethodError volt_to_pressure_and_phase(1e6, :Onda_HGL0200_2322, :my_dummy)
    end

    @testset "enforce_hertz_unit()" begin
        @test enforce_hertz_unit(2.0u"mHz") == 0.002u"Hz"
        @test enforce_hertz_unit(2u"Hz")    == 2u"Hz"
        @test enforce_hertz_unit(2.2u"MHz") == 2.2e6u"Hz"
        @test_throws ErrorException enforce_hertz_unit(2.2u"m")
        @test enforce_hertz_unit(2)     == 2u"Hz"
        @test enforce_hertz_unit(2.2e6) == 2.2e6u"Hz"
    end

    @testset "interp()" begin
        f_axis = collect(0:1:10)
        values = f_axis

        f0_0 = 0.01
        @test interp(f_axis, values, f0_0) == 0.01
        f0_1 = 0.5
        @test interp(f_axis, values, f0_1) == 0.5
        f0_low = -0.1
        @test_throws "f0 must be between $f_axis[1]-$f_axis[end]" interp(f_axis, values, f0_low)
        f0_high = 10.1
        @test_throws "f0 must be between $f_axis[1]-$f_axis[end]" interp(f_axis, values, f0_high)
    end
end