text:
using Plots
#using Pkg; Pkg.add("DistributedArrays")
#using Transducers
using FileIO
using Distributed
using DistributedArrays
using DistributedArrays.SPMD
#@everywhere using SharedArrays
@everywhere using Distributed
@everywhere using DistributedArrays
@everywhere using DistributedArrays.SPMD
using BenchmarkTools

#rmprocs(2:37)
d_closeall()

N_workers = 32
addprocs(N_workers)

println(workers())

@everywhere include_string(Main, $(read("aux_functions.jl", String)), "aux_functions.jl")

order = 1
const L = 5
const N = 80000
#const N = 20000
const J = 50
const κ = 2π/L # for Fourier Transform
dx = L/J
x = [dx*(i-1) for i in 1:J] ;
p = (L, N, J, κ, dx, order)

DS = dzeros((J,nworkers()), workers(), [1,nworkers()]);

pars = zeros(2N)
#pars = dones((N,), workers(), [nworkers()] );
par_dis = load("Initial_Distributions/par_dis_L5_N8_5_theta001.jld2", "par_dis");
reorder_particles!(par_dis,pars)

Dpars = distribute(pars);

if nprocs() > 1
    @benchmark begin
        Sf = zeros(J)
    spmd(get_current_ro_par, Dpars, DS, p)
    for j in 1:J
        for i in 1:nworkers()
            Sf[j] += DS[j,i]
        end
    end
    end
end

if nprocs() == 1
    @benchmark begin
    S = zeros(J)
    get_current_ro!(pars, S, p)
    end
end

#(S - Sf)'*(S-Sf)


text:
@testset "Square root Kalman filter" begin

    @testset "Calculate upper triangular of QR" begin
        A = randn(10, 5)
        R_test = @inferred KalmanFilters.calc_upper_triangular_of_qr!(copy(A))

        Q, R = qr(A)
        @test R_test ≈ R

        qr_zeros = zeros(10)
        qr_space_length = @inferred KalmanFilters.calc_gels_working_size(A, qr_zeros)
        qr_space = zeros(qr_space_length)
        R_res = zeros(5, 5)
        R_test_inplace = @inferred KalmanFilters.calc_upper_triangular_of_qr_inplace!(R_res, copy(A), qr_zeros, qr_space)
        @test R_test_inplace ≈ R
    end

    @testset "Time update" begin

        x = [1., 1.]
        P = [1. 0.; 0. 1.]
        F = [1. 0.; 0. 2.]
        Q = [1. 0.; 0. 1.]

        tu = @inferred time_update(x, cholesky(P), F, cholesky(Q))
        @test @inferred(get_state(tu)) == [1., 2.]
        @test @inferred(get_covariance(tu)) ≈ [2. 0.; 0. 5.]

        x = 1.
        P = 1.
        F = 1.
        Q = 1.

        tu = @inferred time_update(x, cholesky(P), F, cholesky(Q))
        @test @inferred(get_state(tu)) == 1.
        @test @inferred(get_covariance(tu)) ≈ [2.]

        x = randn(6)
        PL = randn(6,6)
        P = PL'PL
        P_chol = cholesky(P)
        QL = randn(6,6)
        Q = QL'QL
        Q_chol = cholesky(Q)
        F = randn(6,6)

        tu = @inferred time_update(x, P, F, Q)
        tu_chol = @inferred time_update(x, P_chol, F, Q_chol)
        @test @inferred(get_covariance(tu_chol)) ≈ get_covariance(tu)
        @test @inferred(get_state(tu_chol)) ≈ get_state(tu)

        tu_interm = @inferred SRKFTUIntermediate(6)
        tu_chol_inplace = @inferred time_update!(tu_interm, x, P_chol, F, Q_chol)
        @test @inferred(get_covariance(tu_chol_inplace)) ≈ get_covariance(tu)
        @test @inferred(get_state(tu_chol_inplace)) ≈ get_state(tu)
    end

    @testset "Measurement update" begin

        y = [1., 1.]
        x = [1., 1.]
        P = [1. 0.; 0. 1.]
        H = [1. 0.; 0. 1.]
        R = [1. 0.; 0. 1.]

        mu = @inferred measurement_update(x, cholesky(P), y, H, cholesky(R))
        @test @inferred(get_state(mu)) == [1., 1.]
        @test @inferred(get_covariance(mu)) ≈ [0.5 0.; 0. 0.5]
        @test @inferred(get_innovation(mu)) == [0.0, 0.0]
        @test @inferred(get_innovation_covariance(mu)) ≈ [2.0 0.0; 0.0 2.0]
        @test @inferred(get_kalman_gain(mu)) ≈ [0.5 0.0; 0.0 0.5]

        x = randn(6)
        PL = randn(6,6)
        P = PL'PL
        P_chol = cholesky(P)
        RL = randn(3,3)
        R = RL'RL
        y = randn(3)
        R_chol = cholesky(R)
        H = randn(3,6)

        mu = @inferred measurement_update(x, P, y, H, R)
        mu_chol = @inferred measurement_update(x, P_chol, y, H, R_chol)
        @test @inferred(get_covariance(mu_chol)) ≈ get_covariance(mu)
        @test @inferred(get_state(mu_chol)) ≈ get_state(mu)

        mu_interm = SRKFMUIntermediate(6, 3)
        mu_chol_inplace = @inferred measurement_update!(mu_interm, x, P_chol, y, H, R_chol)
        @test @inferred(get_covariance(mu)) ≈ get_covariance(mu_chol_inplace)
        @test @inferred(get_state(mu)) ≈ get_state(mu_chol_inplace)
    end
end

ext: jl


text:
"""
    pointcloud

This plot uses matrix and vector of vector both.
"""
@userplot PointCloud
@recipe function f(pc::PointCloud)
    P = pc.args[1]

    x = []
    y = []
    if typeof(P) <: Matrix
        x = P[1,:]
        y = P[2,:]
    elseif typeof(P) <: Vector
        for p in P
            push!(x, p[1])
            push!(y, p[2])
        end
    end
    @series begin
        # label --> :none
        seriestype := :scatter
        x, y
    end
end


text:
__precompile__()
module SenseHat

export led_matrix, RGB565, led_clear,
    show_char, show_message,
    Stick, StickEvent, readstick, sticktask,
    pressure, temperature, humidity,
    @colorant_str

include("ioctl.jl")

include("led.jl")
using .LED
include("led_extra.jl")

include("stick.jl")
using .Stick

include("sensors.jl")
using .Sensors

include("text.jl")
using .Text

end # module


text:
# construct graph
G = DiGraph(3)
add_edge!(G,1,2)
add_edge!(G,2,3)

# source and sink terminals
source, target = 1, 3

# default capacity
capacity_matrix = LightGraphs.DefaultCapacity(G)

# state variables
flow_matrix = zeros(Int, 3, 3)

TREE = zeros(Int, 3)
TREE[source] = 1
TREE[target] = 2

PARENT = zeros(Int, 3)

A = [source,target]

residual_graph = LightGraphs.residual(G)

path = LightGraphs.find_path!(residual_graph, source, target, flow_matrix, capacity_matrix, PARENT, TREE, A)

@test path == [1,2,3]


text:
module RouteguideClients
using gRPCClient

include("routeguide.jl")
using .routeguide

import Base: show

# begin service: routeguide.RouteGuide

export RouteGuideBlockingClient, RouteGuideClient

struct RouteGuideBlockingClient
    controller::gRPCController
    channel::gRPCChannel
    stub::RouteGuideBlockingStub

    function RouteGuideBlockingClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = RouteGuideBlockingStub(channel)
        new(controller, channel, stub)
    end
end

struct RouteGuideClient
    controller::gRPCController
    channel::gRPCChannel
    stub::RouteGuideStub

    function RouteGuideClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = RouteGuideStub(channel)
        new(controller, channel, stub)
    end
end

show(io::IO, client::RouteGuideBlockingClient) = print(io, "RouteGuideBlockingClient(", client.channel.baseurl, ")")
show(io::IO, client::RouteGuideClient) = print(io, "RouteGuideClient(", client.channel.baseurl, ")")

import .routeguide: GetFeature
"""
    GetFeature

- input: routeguide.Point
- output: routeguide.Feature
"""
GetFeature(client::RouteGuideBlockingClient, inp::routeguide.Point) = GetFeature(client.stub, client.controller, inp)
GetFeature(client::RouteGuideClient, inp::routeguide.Point, done::Function) = GetFeature(client.stub, client.controller, inp, done)

import .routeguide: ListFeatures
"""
    ListFeatures

- input: routeguide.Rectangle
- output: Channel{routeguide.Feature}
"""
ListFeatures(client::RouteGuideBlockingClient, inp::routeguide.Rectangle) = ListFeatures(client.stub, client.controller, inp)
ListFeatures(client::RouteGuideClient, inp::routeguide.Rectangle, done::Function) = ListFeatures(client.stub, client.controller, inp, done)

import .routeguide: RecordRoute
"""
    RecordRoute

- input: Channel{routeguide.Point}
- output: routeguide.RouteSummary
"""
RecordRoute(client::RouteGuideBlockingClient, inp::Channel{routeguide.Point}) = RecordRoute(client.stub, client.controller, inp)
RecordRoute(client::RouteGuideClient, inp::Channel{routeguide.Point}, done::Function) = RecordRoute(client.stub, client.controller, inp, done)

import .routeguide: RouteChat
"""
    RouteChat

- input: Channel{routeguide.RouteNote}
- output: Channel{routeguide.RouteNote}
"""
RouteChat(client::RouteGuideBlockingClient, inp::Channel{routeguide.RouteNote}) = RouteChat(client.stub, client.controller, inp)
RouteChat(client::RouteGuideClient, inp::Channel{routeguide.RouteNote}, done::Function) = RouteChat(client.stub, client.controller, inp, done)

# end service: routeguide.RouteGuide

end # module RouteguideClients


text:
import Base: show, getindex, setindex!, getproperty, setproperty!
import MATLAB: mxarray, put_variable
using MATLAB

# Actual ints are rare in MATLAB: integers are usually stored as doubles.
# Sometimes you get into trouble if the datatype of an integer value
# is actually an int.
mxarray(a :: Int) = mxarray(Float64(a))

# equivalent of a function handle in MATLAB
mutable struct MatFcn
	funcName :: String
end

show(io :: IO, fcn :: MatFcn) = print(io, "@"*fcn.funcName)

# make an options type whose syntax is very similar to that of MATLAB structs
struct Options
	dict :: Dict{String, Any}
	Options() = new(Dict([]))
end

function show(io :: IO, o :: Options)
	print(io, "Options consisting of ")
	show(io, MIME"text/plain"(), o.dict)
end

getindex(options :: Options, key...) = getindex(options.dict, key...)
setindex!(options :: Options, value, key...) = setindex!(options.dict, value, key...)
getproperty(options :: Options, name :: Symbol) = name == :dict ? getfield(options, name) : getindex(options.dict, string(name))
setproperty!(options :: Options, name :: Symbol, x) = setindex!(options.dict, x, string(name))

# When we put a MatFcn into MATLAB, make it serve as a function handle
function convertMatFcn(sess :: MSession, name :: Symbol, fcn :: MatFcn)
	eval_string(sess, string(name, " = str2func('", fcn.funcName, "');"))
end

# Put the MatFcn into MATLAB as a function handle
function put_variable(sess :: MSession, name :: Symbol, fcn :: MatFcn)
	put_variable(sess, name, fcn.funcName)
	convertMatFcn(sess, name, fcn)
end

# Put the Options into MATLAB where all the subfields of type Options and
# MatFcn are converted correctly
function put_variable(sess :: MSession, name :: Symbol, options :: Options)
	put_variable(sess, name, options.dict)
	for (option, value) in options.dict
		if typeof(value) <: MatFcn
			convertMatFcn(sess, Symbol(string(string(name), ".", option)), value)
		elseif typeof(value) <: Options
			# save recursively
			tmpname = "jl_temp_var"
			put_variable(sess, Symbol(tmpname), value)
			eval_string(sess, string(name, ".", option, " = ", tmpname, ";"))
			eval_string(sess, string("clear ", tmpname))
		end
	end
end


text:
abstract type Timestep end

struct FixedTimestep{T} <: Timestep
    dt::T
end

struct RandomTimestep{T} <: Timestep
    distribution::T
end

get_step(timestep::FixedTimestep) = timestep.dt
get_step(timestep::RandomTimestep) = rand(timestep.distribution)


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ ba447bcc-0094-11eb-068b-6b5e045fd602
begin
	using Pkg
	Pkg.activate(".")
	using RvSpectMLBase
	import EchelleInstruments
	import EchelleInstruments: EXPRES
	using DataFrames, Query
	using Statistics
	using Plots
	using PlutoUI
end

# ╔═╡ 78bc57ea-0140-11eb-2f53-15d0cb1e8b7b
md"# Interactive plots of 2 spectra"

# ╔═╡ 19fc60d6-009a-11eb-3d9c-0f9152ece864
md"""Divide by blaze: $(@bind div_blaze CheckBox())   Log y-axis? $(@bind log_y CheckBox())  Show tellurics? $(@bind plt_tellurics CheckBox())  Show continuum? $(@bind plt_continuum CheckBox())    $(@bind go_reset_range Button("Reset pixel range")) """

# ╔═╡ 53b7cab4-04e8-11eb-33c2-65008a5666a4
md"### Plotting functions"

# ╔═╡ 9d4a5640-00a3-11eb-11f7-5d5e9c00fc61
md"### Boring code below"

# ╔═╡ ce23b456-00a5-11eb-0d07-2b15dfff6257
md"### Find files to read data from"

# ╔═╡ c5b8c394-013d-11eb-2353-7da0efded8f9
begin
	expres_data_path = joinpath(homedir(),"Data/EXPRES/inputs")
	target_subdir = "101501"
 	fits_target_str = "101501"
 	output_dir = "examples/output/"
    df_files = EchelleInstruments.make_manifest(expres_data_path, target_subdir, EXPRES );
end;

# ╔═╡ c364ee3e-0095-11eb-0faa-656ede64325f
begin
	linelist_for_ccf_filename = "G8.espresso.mas"
	ccf_mid_velocity = -5e3
	max_spectra_to_use = 100
	df_files_use = df_files |>
	    @filter( _.target == fits_target_str ) |>
	    @take(max_spectra_to_use) |>
	    DataFrame;
end;

# ╔═╡ b00d8782-04f1-11eb-2c87-9d93fc5ec593
begin
	using FITSIO
	f = FITS(df_files_use.Filename[1])
	read_header(f[3])
end

# ╔═╡ 99e1c70c-00a5-11eb-3cc9-fb8934931026
md"### Read data files & extract each order"

# ╔═╡ 01c2dfa6-0096-11eb-1361-af95574e6f95
all_spectra = map(EXPRES.read_data,eachrow(df_files_use));

# ╔═╡ 1a16d0ec-0095-11eb-3bac-0ddaaf387169
function extract_orders(data::AbstractVector{ST}; verbose::Bool = false, orders_to_use = RvSpectMLBase.orders_to_use_default(first(data).inst) ) where { ST<:AbstractSpectra }
      order_list_timeseries = make_order_list_timeseries(data, orders_to_use=orders_to_use)
      order_list_timeseries = filter_bad_chunks(order_list_timeseries,verbose=true)
end;

# ╔═╡ fdcf688a-0091-11eb-3fb4-9333baf67eab
order_list_timeseries = extract_orders(all_spectra );

# ╔═╡ cfabe62c-0098-11eb-3b6d-ed7ab2f9f0f1
md"Select the observation indices to plot: $(@bind obs_idx Slider(1:length(order_list_timeseries); default=1, show_value=true))   $(@bind obs_jdx Slider(1:length(order_list_timeseries); default=2, show_value=true))"

# ╔═╡ a6031a34-0098-11eb-38aa-891004881de9
md"Select a chunk ID: $(@bind chunk_idx Slider(1:num_chunks(order_list_timeseries); default=floor(Int,num_chunks(order_list_timeseries)//2), show_value=true))   Order = $(order_list_timeseries[obs_idx].order[chunk_idx])"

# ╔═╡ 4ab1c184-009b-11eb-3438-35502bf52add
begin
	go_reset_range
	pixel_range_min = 1
	pixel_range_max = maximum(map(chid->length(order_list_timeseries[1][chid].λ),1:num_chunks(order_list_timeseries)))
	max_width_pixels = floor(Int,(pixel_range_max)//2)
	md"Central pixel: $(@bind mid_pixel Slider(pixel_range_min:pixel_range_max; default=floor(Int,(pixel_range_min+pixel_range_max)//2)))  and
	Half-width in pixels: $(@bind halfwidth_pixels Slider(1:max_width_pixels; default=max_width_pixels))"
end

# ╔═╡ 066347f4-04d9-11eb-3a82-4f986133e7ac
get_meta(key::Symbol, obsid::Integer; digits::Integer = 3) = floor(parse(Float64,order_list_timeseries.metadata[obsid][key])*10^3)/10^3;

# ╔═╡ f226b23a-04d1-11eb-02ed-8d71509a5731
begin
	print_ti = round((order_list_timeseries.times[obs_idx]-order_list_timeseries.times[1])*100)/100
	print_tj = round((order_list_timeseries.times[obs_jdx]-order_list_timeseries.times[1])*100)/100
	pwv1 = get_meta(:PWV,obs_idx)
	pwv2 = get_meta(:PWV,obs_jdx)
	pwv2 = get_meta(:PWV,obs_jdx)
	air1 = get_meta(:AIRMASS,obs_idx)
	air2 = get_meta(:AIRMASS,obs_jdx)
	sun1 = get_meta(:SUNDIST,obs_idx)
	sun2 = get_meta(:SUNDIST,obs_jdx)
	moon1 = get_meta(:MOONDIST,obs_idx)
	moon2 = get_meta(:MOONDIST,obs_jdx)
	bc1 =  NaN # get_meta(:wtd_mdpt_bc,obs_idx)
	bc2 =  NaN # get_meta(:wtd_mdpt_bc,obs_jdx)
end;

# ╔═╡ f416a706-04e7-11eb-3fa1-53dc1bd0401d
string("Times = ", print_ti, " & ", print_tj, "   Δt = ", print_tj-print_ti, " days")

# ╔═╡ 60ccca98-04f1-11eb-0b88-2993a6c13be7
string("BC  = ", bc1, " & ", bc2, "   ΔBC  = ", bc2-bc1)

# ╔═╡ f4d659c0-04e7-11eb-1f1a-05c3603a874d
string("PWV = ", pwv1, " & ", pwv2, "   ΔPWV = ", pwv2-pwv1)

# ╔═╡ e0f528c8-04e7-11eb-135d-c19a28f13d58
string("Airmass = ", air1, " & ", air2, "   ΔAir = ", air2-air1)

# ╔═╡ e60d0470-04e7-11eb-3457-7f4d35f8bffd
string("SunDist  = ", sun1, " & ", sun2, "   ΔD_Sun  = ", sun2-sun1, "    MoonDist = ", moon1, " & ", moon2, "   ΔD_Mood = ", moon2-moon1)

# ╔═╡ c1c3fb94-04cb-11eb-1830-f5a00c592959
md"### Some CSS tweaks to use more of screen"

# ╔═╡ eeb88092-015c-11eb-3989-2de8662a44bd

md"""Display width $(@bind pluto_width_pixels Select([ "600","800","1000","1200","1600","2000"] ))"""

# ╔═╡ 4edfca32-0092-11eb-25d5-dbda1e798503
function compare_two_spectra()
	plt = plot(size=(parse(Int,pluto_width_pixels),700))
	xlabel!("λ (Å)")
	ylabel!("Flux")
	title!("Comparing two spectra")
	# Make copy spectral chunks to plot, since many need to manipulate for plotting
	(pixels1, order1) = order_list_timeseries[obs_idx][chunk_idx].flux.indices
	(pixels2, order2) = order_list_timeseries[obs_jdx][chunk_idx].flux.indices
	flux1 = copy(order_list_timeseries[obs_idx][chunk_idx].flux)
	flux2 = copy(order_list_timeseries[obs_jdx][chunk_idx].flux)
	# Choose which pixels to plot
	min_pixel = min(max(1,mid_pixel-halfwidth_pixels),length(order_list_timeseries[obs_idx][chunk_idx].flux))
	max_pixel = max(min(mid_pixel+halfwidth_pixels, length(order_list_timeseries[obs_idx][chunk_idx].flux)),1)
	# Identify relevant blaze function and continuum
	blaze1 = order_list_timeseries.metadata[obs_idx][:blaze][pixels1,order1]
	blaze2 = order_list_timeseries.metadata[obs_jdx][:blaze][pixels2,order2]
	if plt_continuum
		continuum1 = copy(order_list_timeseries.metadata[obs_idx][:continuum][pixels1,order1])
		continuum2 = copy(order_list_timeseries.metadata[obs_jdx][:continuum][pixels2,order2])
	end
	# Either divide flux by blaze function or multiply continue by blaze
	if div_blaze
		flux1 ./= blaze1
		flux2 ./= blaze2
	else
		if plt_continuum
			continuum1 .*= blaze1
			continuum2 .*= blaze2
		end
	end
	# Calc scaling factor to make spectra line up
	scalefactor = mean(flux1)/mean(flux2)
	flux2 .*= scalefactor
	if plt_continuum
		continuum2 .*= scalefactor
	end
	if log_y
		flux1 = log10.(flux1)
		flux2 = log10.(flux2)
		if plt_continuum
			continuum1 = log10.(continuum1)
			continuum2 = log10.(continuum2)
		end
	end
	# Start plotting stuff
	plot!(plt,order_list_timeseries[obs_idx][chunk_idx].λ[min_pixel:max_pixel],flux1[min_pixel:max_pixel], color=:blue, label=:none) # string(obs_idx))
	plot!(plt,order_list_timeseries[obs_jdx][chunk_idx].λ[min_pixel:max_pixel],flux2[min_pixel:max_pixel], color=:red, label=:none) # string(obs_jdx))
	if plt_tellurics
		tellurics1 = order_list_timeseries.metadata[obs_idx][:tellurics][pixels1,order1]
		tellurics_idx1 = findall(x->x!=one(x),tellurics1[min_pixel:max_pixel]) .+ (min_pixel-1)
		scatter!(plt,order_list_timeseries[obs_idx][chunk_idx].λ[tellurics_idx1],flux1[tellurics_idx1], markersize=1.2, color=:green, label=:none)
		tellurics2 = order_list_timeseries.metadata[obs_idx][:tellurics][pixels2,order2]
		tellurics_idx2 = findall(x->x!=one(x),tellurics2[min_pixel:max_pixel]) .+ (min_pixel-1)
		scatter!(plt,order_list_timeseries[obs_jdx][chunk_idx].λ[tellurics_idx2],flux2[tellurics_idx2], markersize=1.2, color=:green, label=:none)
	end
	if plt_continuum
		plot!(plt,order_list_timeseries[obs_idx][chunk_idx].λ[min_pixel:max_pixel],continuum1[min_pixel:max_pixel], color=:cyan, label=:none)
		plot!(plt,order_list_timeseries[obs_jdx][chunk_idx].λ[min_pixel:max_pixel],continuum2[min_pixel:max_pixel], color=:magenta, label=:none)
	end
	#=
	if plt_mask
		for x in lines_in_mask
			plot([x,x],[ylims(plt)[1], ylims(plt)[2]],color=:black, label=:none)
		end
	end
	=#
	plt
end;

# ╔═╡ b6a6ad72-04d7-11eb-0bd5-81c1ccc85a5c
compare_two_spectra()

# ╔═╡ 9690e8b2-04d8-11eb-1449-a78cf2d8f32f
HTML("<style> main { max-width:" * string(pluto_width_pixels) * "px; } </style> ")

# ╔═╡ d22660c6-015b-11eb-0f6d-8745e4cd4c35
HTML("""<style> input[type="range"] { width:240px; } </style>""")

# ╔═╡ Cell order:
# ╟─78bc57ea-0140-11eb-2f53-15d0cb1e8b7b
# ╟─cfabe62c-0098-11eb-3b6d-ed7ab2f9f0f1
# ╟─a6031a34-0098-11eb-38aa-891004881de9
# ╟─4ab1c184-009b-11eb-3438-35502bf52add
# ╟─19fc60d6-009a-11eb-3d9c-0f9152ece864
# ╟─b6a6ad72-04d7-11eb-0bd5-81c1ccc85a5c
# ╟─f416a706-04e7-11eb-3fa1-53dc1bd0401d
# ╟─60ccca98-04f1-11eb-0b88-2993a6c13be7
# ╟─f4d659c0-04e7-11eb-1f1a-05c3603a874d
# ╟─e0f528c8-04e7-11eb-135d-c19a28f13d58
# ╟─e60d0470-04e7-11eb-3457-7f4d35f8bffd
# ╟─53b7cab4-04e8-11eb-33c2-65008a5666a4
# ╟─4edfca32-0092-11eb-25d5-dbda1e798503
# ╟─066347f4-04d9-11eb-3a82-4f986133e7ac
# ╟─f226b23a-04d1-11eb-02ed-8d71509a5731
# ╟─9d4a5640-00a3-11eb-11f7-5d5e9c00fc61
# ╟─ba447bcc-0094-11eb-068b-6b5e045fd602
# ╟─ce23b456-00a5-11eb-0d07-2b15dfff6257
# ╟─c5b8c394-013d-11eb-2353-7da0efded8f9
# ╟─c364ee3e-0095-11eb-0faa-656ede64325f
# ╟─99e1c70c-00a5-11eb-3cc9-fb8934931026
# ╟─01c2dfa6-0096-11eb-1361-af95574e6f95
# ╟─1a16d0ec-0095-11eb-3bac-0ddaaf387169
# ╠═fdcf688a-0091-11eb-3fb4-9333baf67eab
# ╠═b00d8782-04f1-11eb-2c87-9d93fc5ec593
# ╟─c1c3fb94-04cb-11eb-1830-f5a00c592959
# ╟─eeb88092-015c-11eb-3989-2de8662a44bd
# ╟─9690e8b2-04d8-11eb-1449-a78cf2d8f32f
# ╟─d22660c6-015b-11eb-0f6d-8745e4cd4c35


text:
using Documenter, SearchLight, Validation, PostgreSQLDatabaseAdapter

push!(LOAD_PATH,  "../../src")
push!(LOAD_PATH,  "../../src/adapters")

makedocs()


text:
#=
This file is part of Pseudospectra.jl.

Julia translation
Copyright © 2017 Ralph Smith

Portions from original MATLAB code (via EigTool)
Copyright © 2002-2014 James Burke, Adrian Lewis, Emre Mengi and Michael Overton

SPDX-License-Identifier: BSD-3-Clause
License-Filename: LICENSES/BSD-3-Clause_Eigtool
=#

"""
    psa_abscissa(A,ϵ [,d]) -> α,z

Compute ϵ-pseudospectral abscissa for a dense matrix.

Quadratically convergent two-way method to compute the
ϵ-pseudospectral abscissa `α` of a dense matrix `A`. Also returns a vector
`z` of points where the pseudospectrum reaches the abscissa.
Uses the criss-cross algorithm of Burke et al.

The ϵ-pseudospectral abscissa is
```
   maximum(real(z)) for z s.t. minimum(σ(A-zI)) == ϵ
```

Optional arg:

* `d`: eigenvalues of A, if known in advance

Keyword args:

* `verbosity: 0 for quiet, 1 for noise
"""
function psa_abscissa end

#=
Undocumented keyword args:
* `plotfig`: if a +ve integer, a figure number to use for diagnostic plots (**uses PyPlot**)
* `iterprompt`: whether to pause between steps; for plot visibility

Plotting only works if PyPlot bindings are imported into this module, which we
won't do by default.
=#

function psa_abscissa(A,epsln,eA=[];verbosity=0,plotfig=0,iterprompt=false)
# A = ps_data.input_matrix

# This is the version from EigTool.
#  It is suited to epsilon not too small (otherwise see Overton's web
#  pages, or perhaps Mengi's) and moderate size dense A
#  For large (esp. sparse) A, see Overton's psapsr code.

    if (plotfig > 0)
        if (!isdefined(Main, :PyPlot) || !isdefined(:plot)
            || !(plot === Main.PyPlot.plot))
            @warn("plotting is only implemented for PyPlot")
            plotfig = 0
        end
    end

    n,m = size(A)
    (n==m) || throw(ArgumentError("matrix must be square"))
    (isa(epsln,Real) && (epsln >= 0)) || throw(ArgumentError("ϵ must be >= 0"))

    isempty(eA) && (eA = eigfact(A)[:values])

    if plotfig > 0
        figure(plotfig)
        clf()
        plot(real(eA),imag(eA),"k.")
    end

    smalltol = 1e-10 * max(norm(A),epsln)
    purespa = maximum(real.(eA)) # spectral abscissa

    if epsln == 0 # degenerate case
        idx = sortperm(-real.(eA))
        eA = eA[idx]
        f = purespa
        z = [ew for ew in eA if (real(ew) >= (f-smalltol))]
        return f,z
    else
        if (epsln < 1e-9)
            # TODO: divert to specialized Hamiltonian method
            @warn("epsln too small; expect poor accuracy")
        end
        xold = -Inf
        x,ind = findmax(real.(eA)) # initial iterate
        y = imag(eA[ind])
        (verbosity > 0) && @printf("\npsa_abscissa: x=%22.15f   ",x)
        iter = 0
        no_imageig = false
        ybest = NaN
        E = Matrix(epsln * I, size(A)...)
        Areal = (eltype(A) <: Real)
        imagtol = smalltol # used to detect zero real parts

        while !no_imageig && (x > xold)
            iter += 1
            (iter > 20) && error("psa_abscissa: too many steps")
            yold = y
            # given current x, look for all relevant interseections of vertical
            # line w/ pseudospectrum, process and return pair midpoints.
            # note scalar x, scalar ybest, vector y
            # (updates imagtol)
            y,newitol = pspa_2way_imag(A,E,epsln,x,ybest,iter,imagtol,plotfig,verbosity)
            imagtol = newitol
            ptout = "Computing Pseudospectral Abscissa..(iteration $iter)"
            if iterprompt
                print(" [RET]")
                readline()
            end
            no_imageig = isempty(y)
            if !no_imageig
                (verbosity > 0) && print("psa_abscissa: y = $y  ")

                # given resulting y values, look for rightmost intersection of
                # corresponding horizontal lines w/ pseudospectrum
                # note: scalar x, vector y, scalar ybest (even if tied by conj. pair

                xold = x
                ybestt = ybest
                x,ybest = pspa_2way_real(A,E,y,imagtol,plotfig,xold)
                if verbosity > 0
                    println()
                    print("psa_abscissa: x = $x   ")
                end
                if x < xold
                    x = xold # terminates while loop
                    ybest = ybestt
                    if verbosity > 0
                        println()
                        print("psa_abscissa: could not find bigger x")
                    end
                end # if x < xold
            end # if !no_imageig
        end # while

        if isempty(ybest)
            # TODO: fall back to Hamiltonian scheme
            error("Failed in first iteration")
        end

        f = x
        z = x + ybest*im
        if ybest != 0 # FIXME: abs(ybest) < tol
            return f,[z,conj(z)]
        else
            return f,[z]
        end
    end
end

function pspa_2way_real(A,E,y,imagtol,plotfig,xold)

    n = size(A,1)
    Aisreal = (eltype(A) <: Real)
    xnew = zeros(size(y))
    for j in 1:length(y)
        B2 = A - y[j]*im*I
        M2 = vcat(hcat(1im*B2',E),hcat(-E,1im*B2))
        eM2 = eigvals(M2)
        println("test ",minimum(abs.(real.(eM2)))," $imagtol")
        if minimum(abs.(real.(eM2))) <= imagtol # check for real eigenvalue
            xnew[j] = maximum([imag(ew) for ew in eM2 if (abs(real(ew)) < imagtol)]
                              )
        else
            xnew[j] = -Inf
            error("horizontal search failed: no intersection points found for one of the midpoints")
            # TODO: fall back to alternate version
        end
        if plotfig > 0
            figure(plotfig)
            hold(true)
            plot([xold, xnew[j]],y[j]*ones(2),"m-") # horiz line
            plot([xnew[j]],[y[j]],"b+") # right end point
            if Aisreal
                plot([xold,xnew[j]], -y[j]*ones(2),"m-")
                plot([xnew[j]],[-y[j]],"b+")
            end
        end
    end
    xbest, ind = findmax(xnew) # furthest to right
    ybest = y[ind]
    if plotfig > 0
        plot(xbest,ybest,"b*")
    end
    return xbest, ybest
end

"""
 Search for intersections between vertical line with given x component
 and the pseudospectrum.  Start by looking for imaginary eigenvalues of
 Hamiltonian matrix; if there are none, return the empty vector for `ynew`.
 Otherwise, remove any non-extreme imaginary eigenvalues that don't correspond
 to the smallest singular value.  Then sort the eigenvalues into pairs.
 As do so, ensure that `ywant` (a scalar) is in the eigenvalue list, and if
 not, add an extra pair above and below `ywant` (unless the pair containing
 `ywant` is already a very short interval, indicating near convergence to a
 maximizer).  The omission must be caused by rounding, and if we overlook
 this, the process could terminate to a local minimizer of the
 pseudospectrum instead of a global maximizer (consider minus
 the 5,5 Demmel matrix with `ϵ = 0.01`).
 Finally, return the pair midpoints; there must be at least one.

 Updates `imagtol`
"""
function pspa_2way_imag(A,E,epsln,x,ywant,iter,imagtol,plotfig,verbosity)
    svd_check = true
    Areal = (eltype(A) <: Real)
    n = size(A,1)
    B = A - x*I
    # compute eigenvalues of Hamiltonian matrix M = [-B' E; -E B]
    # where E = ϵ I
    M = vcat(hcat(-B',E),hcat(-E, B))
    eM = eigvals(M)
    minreal = minimum(abs.(real.(eM)))
    if (iter == 1) && (minreal > imagtol)
        imagtol += minreal
    end

    newimagtol = imagtol
    ynew = zeros(0)
    if minreal <= imagtol # check if M has an imaginary eigenvalue
        y = sort(imag.([ew for ew in eM if abs(real(ew)) <= imagtol]))
        ny = length(y)
        if svd_check
            indx2 = [1] # check out non-extreme imaginary parts
            for check in 2:(ny-1)
                j = check
                Ashift = A - (x + im*y[j])*I
                u,s,v = svd(Ashift)
                minval,minind = findmin(abs.(s-epsln))
                if minind == n
                    push!(indx2,j)
                end
            end
            push!(indx2,ny)
            removed = ny - length(indx2)
            if removed > 0
                if verbosity > 0
                    println()
                    print("pspa_2way_imag: singular value test removed $removed eigenvalues")
                end
                y = y[indx2]
            end
        end
        if isodd(length(y))
            # TODO: fallback to better version
            error("odd number of intersection points found by vertical search")
        end
        npairs = floor(Int,length(y) / 2)
        if npairs == 0
            # TODO: fallback to better version
            error("all pairs eliminated by singular value test")
        end

        # organize in pairs and take midpoints
        for j in 1:npairs
            ylow = y[2*j-1]
            yhigh = y[2*j]
            # before taking midpoint, if interval is not very short,
            # check and see if ywant is in this interval, well away from the
            # end points. If so, break pair into two pairs, one above and
            # one below ywant, and take midpoints of both.
            inttol = 0.01 * (yhigh - ylow)
            if (ywant > ylow + inttol) && (ywant < yhigh - inttol)
                push!(ynew,(ylow+ywant)/2)
                push!(ynew,(yhigh+ywant)/2)
            else
                push!(ynew,(ylow+yhigh)/2)
            end
        end
        if plotfig > 0
            figure(plotfig)
            plot(x*ones(2),[maximum(y),minimum(y)],"r-") # vertical line
            plot(x*ones(length(y)),y,"g+") # intersections
            plot(x*ones(length(ynew)),ynew,"bx")
        end

        if Areal
            ynew = [yy for yy in ynew if yy >= 0]
        end
    end
    return ynew, imagtol
end


text:
doc"""
Compute a rigorous bound for a Taylor1Model.
"""
bound(f::Taylor1Model) = bound(f, f.I)

doc"""
Compute a rigorous bound for a `Taylor1Model` `f` evaluated over the interval `I`.
This interval must be a subset of `f.I`.
"""
function bound(f::Taylor1Model, I)
    if !(I ⊆ f.I)
        error("Cannot evaluate `Taylor1Model` on interval $I that is not included in interval of definition, $(f.I)")
    end

    return bound(f.p, f.x0, I) + f.Δ
end

function bound(p::Taylor1, x0, I)

    B = zero(I)
    n = degree(p)

    for i = n:-1:0
        B = B * (I - x0) + bound(p[i])
    end

    return B
end

bound(f::TaylorN) = evaluate(f)  # can replace by better polynomial bounder

# bound(f::Taylor1Model) = bound(f.p, f.x0, f.I)

bound(x::Interval) = x


text:
function parse_kw(ps::ParseState, ::Type{Val{Tokens.MACRO}})
    startbyte = ps.t.startbyte
    start_col = ps.t.startpos[2] + 4
    kw = INSTANCE(ps)
    format_kw(ps)
    if ps.nt.kind == Tokens.IDENTIFIER
        next(ps)
        sig = INSTANCE(ps)
        @catcherror ps startbyte sig = parse_call(ps, sig)
    else
        @catcherror ps startbyte sig = @closer ps block @closer ps ws parse_expression(ps)
    end

    _get_sig_defs!(sig)
    block = EXPR{Block}(EXPR[], 0, Variable[], "")
    @catcherror ps startbyte @default ps parse_block(ps, block, start_col)

    next(ps)
    ret = EXPR{Macro}(EXPR[kw, sig, block, INSTANCE(ps)], ps.nt.startbyte - startbyte, Variable[], "")
    ret.defs =  [Variable(Symbol("@", Expr(_get_fname(sig))), :Macro, ret)]
    return ret
end

"""
    parse_macrocall(ps)

Parses a macro call. Expects to start on the `@`.
"""
function parse_macrocall(ps::ParseState)
    startbyte = ps.t.startbyte
    next(ps)
    mname = EXPR{IDENTIFIER}(EXPR[], ps.nt.startbyte - ps.lt.startbyte, Variable[], string("@", ps.t.val))
    # Handle cases with @ at start of dotted expressions
    if ps.nt.kind == Tokens.DOT && isemptyws(ps.ws)
        while ps.nt.kind == Tokens.DOT
            next(ps)
            op = INSTANCE(ps)
            if ps.nt.kind != Tokens.IDENTIFIER
                return EXPR{ERROR}(EXPR[], 0, Variable[], "Invalid macro name")
            end
            next(ps)
            nextarg = INSTANCE(ps)
            mname = EXPR{BinarySyntaxOpCall}(EXPR[mname, op, Quotenode(nextarg)], mname.span + op.span + nextarg.span, Variable[], "")
        end
    end
    ret = EXPR{MacroCall}(EXPR[mname], 0, Variable[], "")

    if ps.nt.kind == Tokens.COMMA
        ret.span = ps.nt.startbyte - startbyte
        return ret
    end
    if isemptyws(ps.ws) && ps.nt.kind == Tokens.LPAREN
        next(ps)
        push!(ret.args, INSTANCE(ps))
        @catcherror ps startbyte args = @default ps @nocloser ps newline @closer ps paren parse_list(ps, ret.args)
        append!(ret.args, args)
        next(ps)
        push!(ret.args, INSTANCE(ps))
    else
        insquare = ps.closer.insquare
        @default ps while !closer(ps)
            @catcherror ps startbyte a = @closer ps inmacro @closer ps ws @closer ps wsop parse_expression(ps)
            push!(ret.args, a)
            if insquare && ps.nt.kind == Tokens.FOR
                break
            end
        end
    end
    ret.span = ps.nt.startbyte - startbyte
    return ret
end

ismacro(x) = false
ismacro(x::EXPR{LITERAL{Tokens.MACRO}}) = true
ismacro(x::EXPR{Quotenode}) = ismacro(x.args[1])
function ismacro(x::EXPR{BinarySyntaxOpCall})
    if x.args[2] isa OPERATOR{DotOp,Tokens.DOT}
        return ismacro(x.args[2])
    else
        return false
    end
end


text:
#=
 # 121: 1D Poisson with point charge
 ([source code](SOURCE_URL))

Solve a Poisson equation
```math
- \Delta u = 0
```

in $\Omega=(-1,1)$
with a point charge $Q$ at $x=0$.
=#

module Example121_PoissonPointCharge1D

using Printf

using VoronoiFVM
using ExtendableGrids
using GridVisualize

function main(;nref=0,Plotter=nothing, verbose=false, unknown_storage=:sparse, brea=false)

    ## Create grid in (-1,1) refined around 0
    hmax=0.2/2.0^nref
    hmin=0.05/2.0^nref
    X1=geomspace(-1.0,0.0, hmax,hmin)
    X2=geomspace(0.0,1.0, hmin,hmax)
    X=glue(X1,X2)
    grid=VoronoiFVM.Grid(X)

    ## Edit default region numbers:
    ##   additional boundary region 3 at 0.0
    bfacemask!(grid, [0.0],[0.0],3)
    ## Material 1 left of 0
    cellmask!(grid, [-1.0],[0.0],1)
    ## Material 2 right of 0
    cellmask!(grid, [0.0],[1.0],2)

    Q=0.0

    function flux!(f,u0,edge)
        u=unknowns(edge,u0)
        f[1]=u[1,1]-u[1,2]
    end
    function storage!(f,u,node)
        f[1]=u[1]
    end

    ## Define boundary reaction defining charge
    ## Note that the term  is written on  the left hand side, therefore the - sign
    function breaction!(f,u,node)
        if node.region==3
            f[1]=-Q
        end
    end

   ## Create physics
    physics=VoronoiFVM.Physics(
        flux=flux!,
        storage=storage!,
        breaction=breaction!
    )

    ## Create system
    sys=VoronoiFVM.System(grid,physics,unknown_storage=:dense)

    ##  put potential into both regions
    enable_species!(sys,1,[1,2])

    ## Set boundary conditions

    boundary_dirichlet!(sys,1,1,1.0)
    boundary_dirichlet!(sys,1,2,0.0)

    ## Create a solution array
    inival=unknowns(sys)
    U=unknowns(sys)
    inival.=0

    ## Create solver control info
    control=VoronoiFVM.NewtonControl()
    control.verbose=verbose
    if ispyplot(Plotter)
        Plotter.clf()
    end

    ## Solve and plot for several values of charge
    for q in [0.1,0.2,0.4,0.8,1.6]

        if brea
            ## Charge in reaction term
            Q=q
        else
            ## Charge as boundary condition
            sys.boundary_values[1,3]=q
        end
        solve!(U,inival,sys, control=control)

        ## Plot data
        if ispyplot(Plotter)
            Plotter.grid()
            coord=grid[Coordinates]
            Plotter.plot(coord[1,:],U[1,:],label=@sprintf("Q=%.2f",q))
            Plotter.xlabel("x")
            Plotter.ylabel("\$\\phi\$")
            Plotter.legend(loc="upper right")
            Plotter.pause(1.0e-10)
        end
    end
    return sum(U)
end

function test()
    main()≈20.254591679579015
end
end


text:
const PIPELINE_TABS = (; Load, Filter, Process)
const VISUALIZATION_TABS = (; Spreadsheet, Chart, Pipelines)

struct UI
    pipelinetabs::SimpleList
    visualizationtabs::SimpleList
end

function Base.show(io::IO, ui::UI)
    p, v = map(getkey, ui.pipelinetabs), map(getkey, ui.visualizationtabs)
    print(io, "UI with pipelines $(p) and visualizations $(v)")
end

extract_options(sym::Union{Symbol, AbstractString}) = Symbol(sym), NamedTuple()
extract_options((sym, kwargs)::Pair) = Symbol(sym), kwargs

function concatenate(tabs, names, value)
    pipelines = Any[]
    for entry in names
        name, kwargs = extract_options(entry)
        tab = tabs[name](value; kwargs...)
        push!(pipelines, SimpleDict("key" => string(name), "value" => tab))
        value = output(tab)
    end
    return pipelines
end

"""
    UI(table; pipelinetabs=(:Load, :Filter, :Process), visualizationtabs=(:Spreadsheet, :Chart, :Pipelines))

Generate a `UI` with a given table as starting value. `pipelinetabs` denote the list of
pipeline tabs to include in the user interface. Each tab can take one of the following types:
`:Load, :Filter, :Process`.
Repetitions are allowed, for example setting
`tabs=(:Load, :Filter, :Process, :Filter)` would generate a `UI` that
allowes filtering both before and after processing the data.
`visualizationtabs` includes the list of visualization tabs to be included in the UI.
Possible values are `:Spreadsheet`, `Chart` and `:Pipelines`.
"""
function UI(table; pipelinetabs=keys(PIPELINE_TABS), visualizationtabs=keys(VISUALIZATION_TABS))
    obs = Observable(SimpleTable(table))
    pipelines = concatenate(PIPELINE_TABS, pipelinetabs, obs)
    visualizations = Any[]
    for entry in visualizationtabs
        name, kwargs = extract_options(entry)
        visualization = SimpleDict(
            "key" => string(name),
            "value" => VISUALIZATION_TABS[name](getvalue.(pipelines); kwargs...)
        )
        push!(visualizations, visualization)
    end
    return UI(pipelines, visualizations)
end

function jsrender(session::Session, ui::UI)
    evaljs(session, js"document.body.classList.add('bg-gray-100');")
    # manually load all dependencies
    for dep in AllDeps
        JSServe.push!(session, dep)
    end
    pipelinetabs, visualizationtabs = Tabs(ui.pipelinetabs), Tabs(ui.visualizationtabs)
    layout = DOM.div(
            class="grid grid-cols-5 h-full",
            DOM.div(class="col-span-2 pl-8", pipelinetabs),
            DOM.div(class="col-span-3 pl-12 pr-8", visualizationtabs)
        )
    return jsrender(session, layout)
end

"""
    app(table; pipelinetabs=(:Load, :Filter, :Process), visualizationtabs=(:Spreadsheet, :Chart))

Generate a [`UI`](@ref) with a given `table` and lists of pipeline and visualization tabs.
Launch the output as a local app.
"""
function app(table; pipelinetabs=keys(PIPELINE_TABS), visualizationtabs=keys(VISUALIZATION_TABS))
    return App() do
        return DOM.div(AllCSS..., UI(table; pipelinetabs, visualizationtabs))
    end
end

"""
    serve(table;
          pipelinetabs=(:Load, :Filter, :Process),
          visualizationtabs=(:Spreadsheet, :Chart),
          url=Sockets.localhost, port=8081)

Generate a [`UI`](@ref) with a given `table` and lists of pipeline and visualization tabs.
Serve the output at the given `url` and `port`.
Return a `JSServe.Server` object.
"""
function serve(table; url=Sockets.localhost, port=8081, verbose=true,
               pipelinetabs=keys(PIPELINE_TABS), visualizationtabs=keys(VISUALIZATION_TABS))
    return Server(app(table; pipelinetabs, visualizationtabs), string(url), port; verbose)
end


text:
#===================================================================================================
Some conventions:
    -All C data structs are wrapped in Julia data structs.
    -The lowest level functions should always be exposed.
    -In C functions where output is passed as an argument, we instead return the output, doing
        standard error checking of the status.

TODO
    I've left a lot of these unfinished, want to see how they typically get used first.
===================================================================================================#

"""
ODBC API Function Definitions
By Jacob Quinn, 2016
In general, the ODBC functions are implemented to mirror the C header files (sql.h,sqlext.h,sqltypes.h,sqlucode.h)
A few liberties are taken in utliizing standard Julia functions and idioms
Format:
  * function name
  * URL reference
  * short function description
  * valid const definitions
  * relevant notes
  * working and tested status
  * function definition code

Contents
 * Macros and Utility Functions
 * Handle Functions
 * Connection Functions
 * Resultset Metadata Functions
 * Query Functions
 * Resultset Retrieval Functions
 * DBMS Meta Functions
 * Error Handling and Diagnostics
"""

# TODO these first few seem complicated
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms712400(v=vs.85).aspx"
function SQLDrivers(env::Ptr{Void},
                    dir,
                    driver_desc::Ptr{SQLWCHAR},
                    desclen,
                    desc_length::Ref{SQLSMALLINT},
                    driver_attr::Ptr{SQLWCHAR},
                    attrlen,
                    attr_length::Ref{SQLSMALLINT})
    @odbc(:SQLDriversW,
                (Ptr{Void}, SQLUSMALLINT, Ptr{SQLWCHAR}, SQLSMALLINT, Ref{SQLSMALLINT}, Ptr{SQLWCHAR},
                 SQLSMALLINT, Ref{SQLSMALLINT}),
                env, dir, driver_desc, desclen, desc_length, driver_attr, attrlen, attr_length)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711004(v=vs.85).aspx"
function SQLDataSources(env::Ptr{Void},
                        dir,
                        dsn_desc::Ptr{SQLWCHAR},
                        desclen,
                        desc_length::Ref{SQLSMALLINT},
                        dsn_attr::Ptr{SQLWCHAR},
                        attrlen,
                        attr_length::Ref{SQLSMALLINT})
    @odbc(:SQLDataSourcesW,
                (Ptr{Void}, SQLUSMALLINT, Ptr{SQLWCHAR}, SQLSMALLINT, Ref{SQLSMALLINT}, Ptr{SQLWCHAR},
                 SQLSMALLINT, Ref{SQLSMALLINT}),
                env, dir, dsn_desc, desclen, desc_length, dsn_attr, attrlen, attr_length)
end

#### Handle Functions ####

function SQLAllocHandle(handletype::SQLSMALLINT, parenthandle::Ptr{Void}, handle::Ref{Ptr{Void}})
    @odbc(:SQLAllocHandle,
                (SQLSMALLINT, Ptr{Void}, Ref{Ptr{Void}}),
                handletype, parenthandle, handle)
end
function SQLAllocHandle(parenthandle::AbstractSQLHandle)
    hdl = Ref{Ptr{Void}}()
    htype = handletype(parenthandle)
    @check parenthandle SQLAllocHandle(handletype(parenthandle), parenthandle.ptr, hdl)
    makehandle(parenthandle, hdl[])
end

"""
# SQLFreeHandle
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms710123(v=vs.85).aspx"

Frees resources associated with a specific environment, connection, statement, or descriptor handle
See SQLAllocHandle for valid handle types
"""
function SQLFreeHandle(handletype::SQLSMALLINT,handle::Ptr{Void})
    @odbc(:SQLFreeHandle,
                (SQLSMALLINT, Ptr{Void}),
                handletype, handle)
end
SQLFreeHandle(hdl::AbstractSQLHandle) = SQLFreeHandle(handletype(hdl), hdl.ptr)

function SQLSetEnvAttr{T<:Union{Int,UInt}}(env_handle::Ptr{Void}, attribute::Int, value::T)
    @odbc(:SQLSetEnvAttr,
                (Ptr{Void}, Int, T, Int), env_handle, attribute, value, 0)
end
function SQLSetEnvAttr{T<:Union{Int,UInt}}(env_handle::SQLEnvironement, attribute::Integer, value::T)
    SQLSetEnvAttr(env_handle.ptr, attribute, value)
end

"""
# SQLGetEnvAttr
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms709276(v=vs.85).aspx"

Returns the current setting of an environment attribute
"""
function SQLGetEnvAttr(env::Ptr{Void},attribute::Int,value::Array{Int,1},bytes_returned::Array{Int,1})
    @odbc(:SQLGetEnvAttr,
                (Ptr{Void}, Int, Ptr{Int}, Int, Ptr{Int}),
                env, attribute, value, 0, bytes_returned)
end
function SQLGetEnvAttr(env::SQLEnvironement, attribute::Int, val::Vector{Int},
                       bytes_returned::Vector{Int})
    SQLGetEnvAttr(env.ptr, attribute, val, bytes_returned)
end

"""
# SQLSetConnectAttr
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713605(v=vs.85).aspx"

Sets attributes that govern aspects of connections.
length of string or binary stream
"""
function SQLSetConnectAttr(dbc::Ptr{Void},attribute::Int,value::UInt,value_length::Int)
    @odbc(:SQLSetConnectAttrW,
                (Ptr{Void},Int,UInt,Int),
                dbc,attribute,value,value_length)
end
function SQLSetConnectAttr(dbc::Ptr{Void},attribute::Int,value::Array{Int},value_length::Int)
    @odbc(:SQLSetConnectAttrW,
                (Ptr{Void},Int,Ptr{Int},Int),
                dbc,attribute,value,value_length)
end
function SQLSetConnectAttr(dbc::SQLConnection, attribute::Int, value::Array{Int}, value_length::Int)
    SQLSetConnectAttr(dbc.ptr, attribute, value, value_length)
end

"""
# SQLGetConnectAttr
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms710297(v=vs.85).aspx"

Returns the current setting of a connection attribute.
"""
function SQLGetConnectAttr{T,N}(dbc::Ptr{Void}, attribute::Int, value::Array{T,N},
                                bytes_returned::Vector{Int})
    @odbc(:SQLGetConnectAttrW,
                (Ptr{Void},Int,Ptr{T},Int,Ptr{Int}),
                dbc,attribute,value,sizeof(T)*N,bytes_returned)
end
function SQLGetConnectAttr{T,N}(dbc::SQLConnection, attribute::Int, val::Array{T,N},
                                bytes_returned::Vector{Int})
    SQLGetConnectAttr(dbc.ptr, attribute, val, bytes_returned)
end

"""
# SQLSetStmtAttr
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms712631(v=vs.85).aspx"

Sets attributes related to a statement.
**this sets the rowset size for ExtendedFetch and FetchScroll**
See SQLSetConnectAttr; SQL_IS_POINTER, SQL_IS_INTEGER, SQL_IS_UINTEGER, SQL_NTS
"""
function SQLSetStmtAttr(stmt::Ptr{Void}, attribute::Int, value::Ref{SQLLEN}, value_length::Int)
    @odbc(:SQLSetStmtAttrW,
                (Ptr{Void},SQLINTEGER,Ref{SQLLEN},SQLINTEGER),
                stmt,attribute,value,value_length)
end
function SQLSetStmtAttr(stmt::Ptr{Void}, attribute::Int, value::Int, value_length::Int)
    @odbc(:SQLSetStmtAttrW,
                (Ptr{Void},SQLINTEGER,SQLULEN,SQLINTEGER),
                stmt,attribute,value,value_length)
end
function SQLSetStmtAttr(stmt::SQLStatement, attribute::Int, value::Union{Ref{SQLLEN},Int},
                        value_length::Int)
    SQLSetStmtAttr(stmt.ptr, attribute, value, value_length)
end

"""
# SQLGetStmtAttr

http://msdn.microsoft.com/en-us/library/windows/desktop/ms715438(v=vs.85).aspx
"""
function SQLGetStmtAttr{T,N}(stmt::Ptr{Void}, attribute::Int, value::Array{T,N},
                             bytes_returned::Vector{Int})
    @odbc(:SQLGetStmtAttrW,
                (Ptr{Void},Int,Ptr{T},Int,Ptr{Int}),
                stmt,attribute,value,sizeof(T)*N,bytes_returned)
end
function SQLGetStmtAttr{T,N}(stmt::SQLStatement, attribute::Int, value::Array{T,N},
                             bytes_returned::Vector{Int})
    SQLGetStmtAttr(stmt.ptr, attribute, value, bytes_returned)
end

"""
# SQLFreeStmt
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms709284(v=vs.85).aspx"

Sops processing associated with a specific statement,
closes any open cursors associated with the statement,
discards pending results, or, optionally,
frees all resources associated with the statement handle.
"""
function SQLFreeStmt(stmt::Ptr{Void},param::UInt16)
    @odbc(:SQLFreeStmt,
                (Ptr{Void},UInt16),
                stmt, param)
end
SQLFreeStmt(stmt::SQLStatement, param::UInt16) = SQLFreeStmt(stmt, param)

"""
# SQLSetDescField
http://msdn.microsoft.com/en-us/library/windows/desktop/ms713560(v=vs.85).aspx
"""
function SQLSetDescField{T,N}(desc::Ptr{Void}, i::Int16, field_id::Int16, value::Array{T,N},
                              value_length::Vector{Int})
    @odbc(:SQLSetDescFieldW,
                (Ptr{Void},Int16,Int16,Ptr{T},Int),
                desc,field_id,value,value_length)
end
function SQLSetDescField{T,N}(desc::SQLDescriptor, i::Int16, field_id::Int16, value::Array{T,N},
                              value_length::Vector{Int})
    SQLSetDescField(desc.ptr, i, field_id, value, value_length)
end

"""
# SQLGetDescField

http://msdn.microsoft.com/en-us/library/windows/desktop/ms716370(v=vs.85).aspx
"""
function SQLGetDescField{T,N}(desc::Ptr{Void}, i::Int16, attribute::Int16, value::Array{T,N},
                              bytes_returned::Array{Int,1})
    @odbc(:SQLGetDescFieldW,
                (Ptr{Void},Int16,Int16,Ptr{T},Int,Ptr{Int}),
                desc,attribute,value,sizeof(T)*N,bytes_returned)
end
function SQLGetDescField{T,N}(desc::SQLDescriptor, i::Int16, attribute::Int16, value::Array{T,N},
                              bytes_returned::Vector{Int})
    SQLGetDescField(desc.ptr, i, attribute, value, bytes_returned)
end

"""
# SQLGetDescRec

http://msdn.microsoft.com/en-us/library/windows/desktop/ms710921(v=vs.85).aspx
"""
function SQLGetDescRec(desc::Ptr{Void}, i::Int16, name::Vector{UInt8},
                       name_length::Vector{Int16}, type_ptr::Vector{Int16},
                       subtype_ptr::Vector{Int16}, length_ptr::Vector{Int},
                       precision_ptr::Vector{Int16}, scale_ptr::Vector{Int16},
                       nullable_ptr::Vector{Int16})
    @odbc(:SQLGetDescRecW,
                (Ptr{Void},Int16,Ptr{UInt8},Int16,Ptr{Int16},Ptr{Int16},Ptr{Int16},
                 Ptr{Int},Ptr{Int16},Ptr{Int16},Ptr{Int16}),
                desc,i,name,length(name),name_length,type_ptr,subtype_ptr,length_ptr,
                precision_ptr,scale_ptr,nullable_ptr)
end
function SQLGetDescRec(desc::SQLDescriptor, i::Int16, name::Vector{UInt8},
                       name_length::Vector{Int16}, type_ptr::Vector{Int16},
                       subtype_ptr::Vector{Int16}, length_ptr::Vector{Int},
                       precision_ptr::Vector{Int16}, scale_ptr::Vector{Int16},
                       nullable_ptr::Vector{Int16})
    SQLGetDescRec(desc.ptr, i, name, name_length, type_ptr, subtype_ptr, length_ptr,
                  precision_ptr, scale_ptr, nullable_ptr)
end

"""
# SQLCopyDesc

http://msdn.microsoft.com/en-us/library/windows/desktop/ms715378(v=vs.85).aspx
"""
function SQLCopyDesc(source_desc::Ptr{Void},dest_desc::Ptr{Void})
    @odbc(:SQLCopyDesc,
                (Ptr{Void},Ptr{Void}),
                source_desc,dest_desc)
end
function SQLCopyDesc(source_desc::SQLDescriptor, dest_desc::SQLDescriptor)
    SQLCopyDesc(source_desc.ptr, dest_desc.ptr)
end

### Connection Functions ###

"""
# SQLConnect
http://msdn.microsoft.com/en-us/library/windows/desktop/ms711810(v=vs.85).aspx

establishes connections to a driver and a data source
"""
function SQLConnect(dbc::Ptr{Void}, servername::AbstractString, username::AbstractString,
                    password::AbstractString)
    @odbc(:SQLConnectW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                dbc,transcode(SQLWCHAR,servername),length(transcode(SQLWCHAR,servername)),
                transcode(SQLWCHAR,username),length(transcode(SQLWCHAR,username)),
                transcode(SQLWCHAR,password),length(transcode(SQLWCHAR,password)))
end
function SQLConnect(dbc::SQLConnection, servername::AbstractString, username::AbstractString,
                    password::AbstractString)
    SQLConnect(dbc.ptr, servername, username, password)
end

"""
# SQLDriverConnect

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms715433(v=vs.85).aspx"
"""
function SQLDriverConnect(dbc::Ptr{Void}, window_handle::Ptr{Void}, conn_string,
                          out_conn::Ptr{SQLWCHAR}, out_len, out_buff::Ref{Int16},
                          driver_prompt)
    @odbc(:SQLDriverConnectW,
                (Ptr{Void},Ptr{Void},Ptr{SQLWCHAR},SQLSMALLINT,Ptr{SQLWCHAR},SQLSMALLINT,
                 Ptr{SQLSMALLINT},SQLUSMALLINT),
                dbc,window_handle,transcode(SQLWCHAR,conn_string),
                length(transcode(SQLWCHAR,conn_string)),out_conn,out_len,out_buff,driver_prompt)
end
# TODO continue from here!

#SQLBrowseConnect
 "http://msdn.microsoft.com/en-us/library/windows/desktop/ms714565(v=vs.85).aspx"
 #Description: supports an iterative method of discovering and enumerating the attributes and attribute values required to connect to a data source
 #Status:
function SQLBrowseConnect(dbc::Ptr{Void},instring::AbstractString,outstring::Array{SQLWCHAR,1},indicator::Array{Int16,1})
    @odbc(:SQLBrowseConnectW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{Int16}),
                dbc,transcode(SQLWCHAR,instring),length(transcode(SQLWCHAR,instring)),transcode(SQLWCHAR,outstring),length(transcode(SQLWCHAR,outstring)),indicator)
end
#SQLDisconnect
 "http://msdn.microsoft.com/en-us/library/windows/desktop/ms713946(v=vs.85).aspx"
 #Description: closes the connection associated with a specific connection handle
 #Status:
function SQLDisconnect(dbc::Ptr{Void})
    @odbc(:SQLDisconnect,
                (Ptr{Void},),
                dbc)
end
#SQLGetFunctions
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms709291(v=vs.85).aspx"
#Descriptions:
#Valid functionid

#supported will be SQL_TRUE or SQL_FALSE
#Status:
function SQLGetFunctions(dbc::Ptr{Void},functionid::UInt16,supported::Array{UInt16,1})
    @odbc(:SQLGetFunctions,
                (Ptr{Void},UInt16,Ptr{UInt16}),
                dbc,functionid,supported)
end

#SQLGetInfo
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711681(v=vs.85).aspx"
#Description:
#Status:
function SQLGetInfo{T,N}(dbc::Ptr{Void},attribute::Int,value::Array{T,N},bytes_returned::Array{Int,1})
    @odbc(:SQLGetInfoW,
                (Ptr{Void},Int,Ptr{T},Int,Ptr{Int}),
                dbc,attribute,value,sizeof(T)*N,bytes_returned)
end

#### Query Functions ####
#SQLNativeSql
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms714575(v=vs.85).aspx"
#Description: returns the SQL string as modified by the driver
#Status:
function SQLNativeSql(dbc::Ptr{Void}, query_string::AbstractString,
                      output_string::Array{SQLWCHAR,1}, length_ind::Array{Int,1})
    @odbc(:SQLNativeSql,
                (Ptr{Void},Ptr{SQLWCHAR},Int,Ptr{SQLWCHAR},Int,Ptr{Int}),
                dbc,transcode(SQLWCHAR,query_string),length(transcode(SQLWCHAR,query_string)),
                output_string,length(output_string),length_ind)
end

#SQLGetTypeInfo
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms714632(v=vs.85).aspx"
#Description:
#valid sqltype
#const SQL_ALL_TYPES =
#Status:
function SQLGetTypeInfo(stmt::Ptr{Void},sqltype::Int16)
    @odbc(:SQLGetTypeInfo,
                (Ptr{Void},Int16),
                stmt,sqltype)
end
SQLGetTypeInfo(stmt::SQLStatement, sqltype::Int16) = SQLGetTypeInfo(stmt.ptr, sqltype)

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713824(v=vs.85).aspx"
function SQLPutData{T}(stmt::Ptr{Void},data::Array{T},data_length::Int)
    @odbc(:SQLPutData,
                (Ptr{Void},Ptr{T},Int),
                stmt,data,data_length)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms710926(v=vs.85).aspx"
function SQLPrepare(stmt::Ptr{Void},query_string::AbstractString)
    @odbc(:SQLPrepareW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,query_string),length(transcode(SQLWCHAR,query_string)))
end
SQLPrepare(stmt::SQLStatement, query_string::AbstractString) = SQLPrepare(stmt.ptr, query_string)

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713584(v=vs.85).aspx"
function SQLExecute(stmt::Ptr{Void})
    @odbc(:SQLExecute,
                (Ptr{Void},),
                stmt)
end
SQLExecute(stmt::SQLStatement) = SQLExecute(stmt.ptr)

#SQLExecDirect
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713611(v=vs.85).aspx"
#Description: executes a preparable statement
#Status:
function SQLExecDirect(stmt::Ptr{Void},query::AbstractString)
    @odbc(:SQLExecDirectW,
                (Ptr{Void},Ptr{SQLWCHAR},Int),
                stmt,transcode(SQLWCHAR,query),length(transcode(SQLWCHAR,query)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms714112(v=vs.85).aspx"
function SQLCancel(stmt::Ptr{Void})
    @odbc(:SQLCancel,
                (Ptr{Void},),
                stmt)
end

#### Resultset Metadata Functions ####
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms715393(v=vs.85).aspx"
function SQLNumResultCols(stmt::Ptr{Void},cols::Ref{Int16})
    @odbc(:SQLNumResultCols,
                (Ptr{Void},Ref{Int16}),
                stmt, cols)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711835(v=vs.85).aspx"
function SQLRowCount(stmt::Ptr{Void},rows::Ref{Int})
    @odbc(:SQLRowCount,
                (Ptr{Void},Ref{Int}),
                stmt, rows)
end

# "http://msdn.microsoft.com/en-us/library/windows/desktop/ms713558(v=vs.85).aspx"
# function SQLColAttribute(stmt::Ptr{Void},x::Int,)
# = @odbc(:SQLColAttributeW,
#                 (Ptr{Void},UInt16,UInt16,Ptr,Int16,Ptr{Int16},Ptr{Int}),
#                 stmt,x,)
# end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms716289(v=vs.85).aspx"
function SQLDescribeCol(stmt,x,nm::Ptr{SQLWCHAR},nmlen,len::Ref{Int16},dt::Ref{Int16},cs::Ref{SQLULEN},dd::Ref{Int16},nul::Ref{Int16})
    @odbc(:SQLDescribeColW,
                (Ptr{Void},SQLUSMALLINT,Ptr{SQLWCHAR},SQLSMALLINT,Ref{SQLSMALLINT},Ref{SQLSMALLINT},Ref{SQLULEN},Ref{SQLSMALLINT},Ref{SQLSMALLINT}),
                stmt,x,nm,nmlen,len,dt,cs,dd,nul)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms710188(v=vs.85).aspx"
function SQLDescribeParam(stmt::Ptr{Void},x::Int,sqltype::Array{Int16,1},column_size::Array{Int,1},decimal_digits::Array{Int16,1},nullable::Array{Int16,1})
    @odbc(:SQLDescribeParam,
                (Ptr{Void},UInt16,Ptr{Int16},Ptr{Int},Ptr{Int16},Ptr{Int16}),
                stmt,x,sqltype,column_size,decimal_digits,nullable)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms712366(v=vs.85).aspx"
function SQLParamData(stmt::Ptr{Void},ptr_buffer::Array{Ptr{Void},1})
    @odbc(:SQLParamData,
                (Ptr{Void},Ptr{Void}),
                stmt,ptr_buffer)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms715409(v=vs.85).aspx"
function SQLNumParams(stmt::Ptr{Void},param_count::Array{Int16,1})
    @odbc(:SQLNumParams,
                (Ptr{Void},Ptr{Int16}),
                stmt,param_count)
end

#### Resultset Retrieval Functions ####
#SQLBindParameter
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms710963(v=vs.85).aspx"
#Description:
#Status:
function SQLBindParameter(stmt::Ptr{Void},x::Int,iotype::Int16,ctype::Int16,sqltype::Int16,column_size::Int,decimal_digits::Int,param_value,param_size::Int,len::Ref{SQLLEN})
    @odbc(:SQLBindParameter,
                (Ptr{Void},UInt16,Int16,Int16,Int16,UInt,Int16,Ptr{Void},Int,Ptr{SQLLEN}),
                stmt,x,iotype,ctype,sqltype,column_size,decimal_digits,param_value,param_size,len)
end
SQLSetParam = SQLBindParameter

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711010(v=vs.85).aspx"
function SQLBindCols(stmt::Ptr{Void},x,ctype,mem,jlsize,indicator::Vector{SQLLEN})
    @odbc(:SQLBindCol,
                (Ptr{Void},SQLUSMALLINT,SQLSMALLINT,Ptr{Void},SQLLEN,Ptr{SQLLEN}),
                stmt,x,ctype,mem,jlsize,indicator)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711707(v=vs.85).aspx"
function SQLSetCursorName(stmt::Ptr{Void},cursor::AbstractString)
    @odbc(:SQLSetCursorNameW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,cursor),length(transcode(SQLWCHAR,cursor)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms716209(v=vs.85).aspx"
function SQLGetCursorName(stmt::Ptr{Void},cursor::Array{UInt8,1},cursor_length::Array{Int16,1})
    @odbc(:SQLGetCursorNameW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{Int16}),
                stmt,transcode(SQLWCHAR,cursor),length(transcode(SQLWCHAR,cursor)),cursor_length)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms715441(v=vs.85).aspx"
function SQLGetData(stmt::Ptr{Void},i,ctype,mem,jlsize,indicator::Ref{SQLLEN})
    @odbc(:SQLGetData,
                (Ptr{Void},SQLUSMALLINT,SQLSMALLINT,Ptr{Void},SQLLEN,Ptr{SQLLEN}),
                stmt,i,ctype,mem,jlsize,indicator)
end

#SQLFetchScroll
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms714682(v=vs.85).aspx"
#Description:
#Status:
function SQLFetchScroll(stmt::Ptr{Void},fetch_orientation::Int16,fetch_offset::Int)
    @odbc(:SQLFetchScroll,
                (Ptr{Void},Int16,Int),
                stmt,fetch_orientation,fetch_offset)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713591(v=vs.85).aspx"
function SQLExtendedFetch(stmt::Ptr{Void},fetch_orientation::UInt16,fetch_offset::Int,row_count_ptr::Array{Int,1},row_status_array::Array{Int16,1})
    @odbc(:SQLExtendedFetch,
                (Ptr{Void},UInt16,Int,Ptr{Int},Ptr{Int16}),
                stmt,fetch_orientation,fetch_offset,row_count_ptr,row_status_array)
end

#SQLSetPos
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713507(v=vs.85).aspx"
#Description:
#valid operation
#Status
function SQLSetPos{T}(stmt::Ptr{Void},rownumber::T,operation::UInt16,lock_type::UInt16)
    @odbc(:SQLSetPos,
                (Ptr{Void},T,UInt16,UInt16),
                stmt,rownumber,operation,lock_type)
end #T can be Uint64 or UInt16 it seems

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms714673(v=vs.85).aspx"
function SQLMoreResults(stmt::Ptr{Void})
    @odbc(:SQLMoreResults,
                (Ptr{Void},),
                stmt)
end

#SQLEndTran
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms716544(v=vs.85).aspx"
#Description:
#Status:
function SQLEndTran(handletype::Int16,handle::Ptr{Void},completion_type::Int16)
    @odbc(:SQLEndTran,
                (Int16,Ptr{Void},Int16),
                handletype,handle,completion_type)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms709301(v=vs.85).aspx"
function SQLCloseCursor(stmt::Ptr{Void})
    @odbc(:SQLCloseCursor,
                (Ptr{Void},),
                stmt)
end

#SQLBulkOperations
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms712471(v=vs.85).aspx"
#Description:
#valid operation
const SQL_ADD = UInt16(4) #SQLBulkOperations
const SQL_UPDATE_BY_BOOKMARK = UInt16(5) #SQLBulkOperations
const SQL_DELETE_BY_BOOKMARK = UInt16(6) #SQLBulkOperations
const SQL_FETCH_BY_BOOKMARK = UInt16(7) #SQLBulkOperations
#Status:
function SQLBulkOperations(stmt::Ptr{Void},operation::UInt16)
    @odbc(:SQLBulkOperations,
                (Ptr{Void},UInt16),
                stmt,operation)
end

#### DBMS Meta Functions ####
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711683(v=vs.85).aspx"
function SQLColumns(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,table::AbstractString,column::AbstractString)
    @odbc(:SQLColumnsW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)),transcode(SQLWCHAR,column),length(transcode(SQLWCHAR,column)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms716336(v=vs.85).aspx"
function SQLColumnPrivileges(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,table::AbstractString,column::AbstractString)
    @odbc(:SQLColumnPrivilegesW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)),transcode(SQLWCHAR,column),length(transcode(SQLWCHAR,column)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms709315(v=vs.85).aspx"
function SQLForeignKeys(stmt::Ptr{Void},pkcatalog::AbstractString,pkschema::AbstractString,pktable::AbstractString,fkcatalog::AbstractString,fkschema::AbstractString,fktable::AbstractString)
    @odbc(:SQLForeignKeysW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,pkcatalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,pkschema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,pktable)),transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,fkcatalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,fkschema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,fktable)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711005(v=vs.85).aspx"
function SQLPrimaryKeys(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,table::AbstractString)
    @odbc(:SQLPrimaryKeysW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711701(v=vs.85).aspx"
function SQLProcedureColumns(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,proc::AbstractString,column::AbstractString)
    @odbc(:SQLProcedureColumnsW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),proc,length(proc),transcode(SQLWCHAR,column),length(transcode(SQLWCHAR,column)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms715368(v=vs.85).aspx"
function SQLProcedures(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,proc::AbstractString)
    @odbc(:SQLProceduresW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),proc,length(proc))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711831(v=vs.85).aspx"
function SQLTables(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,table::AbstractString,table_type::AbstractString)
    @odbc(:SQLTablesW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)),table_type,length(table_type))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713565(v=vs.85).aspx"
function SQLTablePrivileges(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,table::AbstractString)
    @odbc(:SQLTablePrivilegesW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)))
end

#SQLStatistics
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711022(v=vs.85).aspx"
#Description:
#Status:
function SQLStatistics(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,table::AbstractString,unique::UInt16,reserved::UInt16)
    @odbc(:SQLStatisticsW,
                (Ptr{Void},Ptr{UInt8},Int16,Ptr{UInt8},Int16,Ptr{UInt8},Int16,UInt16,UInt16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)),unique,reserved)
end

#SQLSpecialColumns
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms714602(v=vs.85).aspx"
#Description:
#Status:
function SQLSpecialColumns(stmt::Ptr{Void},id_type::Int16,catalog::AbstractString,schema::AbstractString,table::AbstractString,scope::Int16,nullable::Int16)
    @odbc(:SQLSpecialColumnsW,
                (Ptr{Void},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Int16,Int16),
                stmt,id_type,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)),scope,nullable)
end

#### Error Handling Functions ####
#TODO: add consts
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms710181(v=vs.85).aspx"
function SQLGetDiagField(handletype::Int16,handle::Ptr{Void},i::Int16,diag_id::Int16,diag_info::Array{SQLWCHAR,1},buffer_length::Int16,diag_length::Array{Int16,1})
    @odbc(:SQLGetDiagFieldW,
                (Int16,Ptr{Void},Int16,Int16,Ptr{SQLWCHAR},Int16,Ptr{Int16}),
                handletype,handle,i,diag_id,transcode(SQLWCHAR,diag_info),buffer_length,transcode(SQLWCHAR,msg_length))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms716256(v=vs.85).aspx"
function SQLGetDiagRec(handletype,handle,i,state::Ptr{SQLWCHAR},native::Ref{SQLINTEGER},error_msg,errlen,msg_length)
    @odbc(:SQLGetDiagRecW,
                (SQLSMALLINT,Ptr{Void},SQLSMALLINT,Ptr{SQLWCHAR},Ref{SQLINTEGER},Ptr{SQLWCHAR},SQLSMALLINT,Ref{SQLSMALLINT}),
                handletype,handle,i,state,native,error_msg,errlen,msg_length)
end


text:
using Test

function reversedims(a::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    permutedims(a, N:-1:1)
end

@testset "test reversedims" begin
    sz=(4,3,2,1)
    dummyX = rand(sz...)
    size(reversedims(dummyX))
    ret=reversedims(dummyX)
    @test reverse(size(dummyX)) == size(ret)
end


text:
"""
    aaft(ts::AbstractArray{T, 1} where T)

Generate a realization of an amplitude adjusted Fourier transform (AAFT)
surrogate ([Theiler et al., 1992](https://www.sciencedirect.com/science/article/pii/016727899290102S)).

# Arguments
 - **`ts`**: the time series for which to generate the surrogate realization.

# References

J. Theiler et al., Physica D *58* (1992) 77-94 (1992).
[https://www.sciencedirect.com/science/article/pii/016727899290102S](https://www.sciencedirect.com/science/article/pii/016727899290102S)

"""
function aaft(ts::AbstractArray{T, 1} where T)
    any(isnan.(ts)) && throw(DomainError(NaN, "The input must not contain NaN values"))
    n = length(ts)

    # Indices that would sort `ts` in ascending order
    ts_sorted = sort(ts)

    # Phase surrogate
    phasesurr = randomphases(ts)

    # Rescale amplitudes according to original time series
    phasesurr[sortperm(phasesurr)] = ts_sorted
    return phasesurr
end


text:
using JuMP, Ipopt
using Gurobi
using DataFrames # for comparison

include(string(pwd(), "/src/struct_network.jl"))
include(string(pwd(), "/src/readnetwork.jl"))
#----------------------------------------------------------------------------------#
solver = "Gurobi"
# solver = "Ipopt"

testsystem = "case_ieee30"

buses, lines, generators, datamat = readnetwork(testsystem)

lineset = 1:length(lines)
busset = 1:length(buses)
genset = 1:length(generators)

B_g = []
for g in genset
   push!(B_g, generators[g].location)
end

B_gn = Array{Array{Int64,1}}(undef, length(buses))
for ii in 1:length(buses)
   B_gn[ii] = Int64[]
end
for g in genset
   push!(B_gn[generators[g].location], g)
end

#---------------------------------------------------------------------------------#
if solver == "Gurobi"
   m = Model(Gurobi.Optimizer)
   #    set_optimizer_attribute(m, "MIPGap", 1e-6)
   # set_optimizer_attribute(m, "Method", 5)
   # set_optimizer_attribute(m, "Presolve", 0) # max 2
   # set_optimizer_attribute(m, "NonConvex", 2)
   #    set_optimizer_attribute(m, "NumericFocus", 3) # max 3
   # set_optimizer_attribute(m, "Threads", 6)
   # set_optimizer_attribute(m, "QCPdual", 1)
elseif solver == "Ipopt"
   m = Model(Ipopt.Optimizer)
   set_optimizer_attribute(m, "max_iter", 1_000_000)
   # set_optimizer_attribute(m, "linear_solver", "ma86")
   # set_optimizer_attribute(m, "print_level", 5) # default = 5 (0~12)
end

@variables(m, begin
   u[b in busset], (start = 1.0) # squared voltage magnitude, |v|^2
   a[l in lineset] >= 0             # squared current magnitude, |I|^2
   fp[l in lineset] # real power flow
   fq[l in lineset] # reactive power flow
   pg[g in genset]  # real power output of generator
   qg[g in genset]  # reactive power output of generator
end)

@constraint(m, PBalance[b in busset],
   sum(fp[l] for l in buses[b].inline)
   -
   sum(fp[l] + lines[l].r * a[l] for l in buses[b].outline)
   +
   sum(pg[g] for g in B_gn[b])
   -
   buses[b].Pd
   -
   buses[b].Gs * u[b]
   ==
   0)

@constraint(m, QBalance[b in busset],
   sum(fq[l] for l in buses[b].inline)
   -
   sum(fq[l] + lines[l].x * a[l] for l in buses[b].outline)
   +
   sum(qg[g] for g in B_gn[b])
   -
   buses[b].Qd
   +
   buses[b].Bs * u[b]
   ==
   0)

@constraint(m, BetweenNodes_reversed[l in lineset],
   u[lines[l].tbus]
   + 2 * (lines[l].r * fp[l] + lines[l].x * fq[l])
   + a[l] * (lines[l].r^2 + lines[l].x^2)
   ==
   u[lines[l].fbus])

@constraint(m, LineCapFW[l in lineset; lines[l].u != 0], (-fp[l])^2 + (-fq[l])^2 <= (lines[l].u)^2) # if lines[l].u == 0, then unlimited line capacity | matpower format
@constraint(m, LineCapBW[l in lineset; lines[l].u != 0], (-fp[l] - a[l] * lines[l].r)^2 + (-fq[l] - a[l] * lines[l].x)^2 <= (lines[l].u)^2)
# @constraint(m, SOCP[l in lineset], ((fp[l])^2 + (fq[l])^2) <= u[lines[l].tbus] * a[l] )
@constraint(m, SOCP[l in lineset], (-2 * fp[l])^2 + (-2 * fq[l])^2 + (u[lines[l].tbus] - a[l])^2 <= (u[lines[l].tbus] + a[l])^2)

# @constraint(m, LineCapFW[l in lineset; lines[l].u != 0], [lines[l].u, -fp[l], -fq[l]] in SecondOrderCone() ) # if lines[l].u == 0, then unlimited line capacity | matpower format
# @constraint(m, LineCapBW[l in lineset; lines[l].u != 0], [lines[l].u, -fp[l] - a[l]*lines[l].r, -fq[l] - a[l]*lines[l].x] in SecondOrderCone())
# @constraint(m, SOCP[l in lineset], [u[lines[l].tbus] + a[l], -2*fp[l], -2*fq[l], u[lines[l].tbus] - a[l]] in SecondOrderCone())

@constraint(m, PGmin[g in genset], pg[g] >= generators[g].Pmin)
@constraint(m, PGmax[g in genset], pg[g] <= generators[g].Pmax)
@constraint(m, QGmin[g in genset], qg[g] >= generators[g].Qmin)
@constraint(m, QGmax[g in genset], qg[g] <= generators[g].Qmax)

@constraint(m, Vmax[b in busset], u[b] <= buses[b].Vmax^2) #upper limit constraint for voltage square
@constraint(m, Vmin[b in busset], u[b] >= buses[b].Vmin^2) #lower limit constraint for voltage square

@expression(m, OCgen[g in genset],
   generators[g].cost[1] * (pg[g] * datamat["baseMVA"])^2
   + generators[g].cost[2] * (pg[g] * datamat["baseMVA"])
   + generators[g].cost[3])
@expression(m, OC, sum(OCgen[g] for g in genset))

@objective(m, Min, OC)

optimize!(m)
@show termination_status(m)
@show objective_value(m)
println()

max_stars_repo_name: JipKim/ACOPF

max_issues_repo_name: JipKim/ACOPF

max_forks_repo_name: JipKim/ACOPF


text:
using Sigma
using DataFrames

## Smokers and Friends: Undirected Graph
## =====================================

# noisyeq(a,b) = flip(ifelse(a == b, 1.0, 0.2))
noisyeq(a,b) = (a == b) | flip(.75)

# Given a list of smokers, and an edge graph
# Create the condition that if you have smoke friends you're more
# likely to be a smoker
function smokers_befriend_smokers(issmoker::RandArray, edges)
  condition = RandVar{Bool}[]
  for i = 1:length(edges), j = 1:length(edges)
    if edges[i,j] == 1
      push!(condition, noisyeq(issmoker[i],issmoker[j]))
    end
  end
  (&)(condition...)
end

# Make condition from observations of who is a smoker and who aint
function make_observes(issmoker, observed_smokers, observed_nonsmokers)
  (&)(vcat([issmoker[i] for i in observed_smokers],
           [!issmoker[i] for i in observed_nonsmokers])...)
end

function smokers(issmoker::RandArray{Bool})
  mrv_condition = smokers_befriend_smokers(issmoker)
  observations = make_observes(issmoker, observed_smokers, observed_nonsmokers)

end

function marginals(issmoker, observed_smokers, observed_nonsmokers)
  marginals = Dict{Int,(RandVar, RandVar)}[]
  for i = 1:nsmokers
    if (i ∉ observed_smokers) && (i ∉ observed_nonsmokers)
      marginals[i] = (issmoker[i], smoker_obs & mrv_condition)
    end
  end
  marginals
end

function test_smokers()
  datapath = joinpath(benchdir, "smokers","problem-7-friends.csv")
  edges = readtable(datapath)
  edges = edges[:,2:end] # remove "Node" column
  nsmokers = length(edges)

  # prior prob of being a smoker
  issmoker = RandArray([flip(.2) for i = 1:length(edges)])

  # observations
  observed_smokers = [1,8,9,15]
  observed_nonsmokers = [4,6]

  mrv_condition = smokers_befriend_smokers(issmoker, edges)
  observations = make_observes(issmoker, observed_smokers, observed_nonsmokers)
  issmoker, mrv_condition & observations
#   marginals()
end

issmoker, condition = test_smokers()
# @show condition.smt
# cond_sample_tlmh(issmoker, condition,1;split=rand_partial_split, solver=dreal3)
# ndims(condition)
# @show cond_sample_mh(issmoker_randvar,mrv_condition,10)

# ## Benchmarks
# # Query1: Compute the marginal probability of smoking for each of the
# # unobserved nodes given the six observed nodes:
# smoker_obs = (&)(vcat([issmoker[i] for i in observed_smokers],[issmoker[i] for i in observed_nonsmokers])...)

# fraction_true(x::Vector{Bool}) = count(identity,x)/length(x)

# marginals = Dict{Int,Float64}[]
# for i = 1:nsmokers
#   if (i ∉ observed_smokers) && (i ∉ observed_nonsmokers)
#     marginals[i] = cond_prob_mh(issmoker[i], smoker_obs & mrv_condition,10)
#   end
# end
# @show marginals


text:
config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5)
optimize!(mock) = MOIU.mock_optimize!(mock, [1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0])
for mock in mocks(optimize!)
    Tests.dsos_univariate_sum_test(mock, config)
end
optimize!(mock) = MOIU.mock_optimize!(mock, [1.0, 1.0, √2, 1.0, 1.0, -√2])
for mock in mocks(optimize!)
    Tests.sdsos_univariate_sum_test(mock, config)
    Tests.sos_univariate_sum_test(mock, config)
end


text:
"An  Intermediate Representation (IR) on steroids."
module OptimizingIR

include("graph.jl")
include("lookup_table.jl")
include("types.jl")
include("passes.jl")
include("build.jl")
include("interpreter.jl")
include("native.jl")
include("print.jl")

end # module


text:
##############################################################################
#
# Text output
#
##############################################################################

function escapedprint(io::IO, x::Any, escapes::AbstractString)
    print(io, x)
end

function escapedprint(io::IO, x::AbstractString, escapes::AbstractString)
    print_escaped(io, x, escapes)
end

function printtable(io::IO,
                    df::AbstractDataFrame;
                    header::Bool = true,
                    separator::Char = ',',
                    quotemark::Char = '"',
                    nastring::AbstractString = "NA")
    n, p = size(df)
    etypes = eltypes(df)
    if header
        cnames = _names(df)
        for j in 1:p
            print(io, quotemark)
            print(io, cnames[j])
            print(io, quotemark)
            if j < p
                print(io, separator)
            else
                print(io, '\n')
            end
        end
    end
    quotestr = string(quotemark)
    for i in 1:n
        for j in 1:p
            if ! (isna(df[j],i))
                if ! (etypes[j] <: Real)
		    print(io, quotemark)
		    escapedprint(io, df[i, j], quotestr)
		    print(io, quotemark)
                else
		    print(io, df[i, j])
                end
            else
		print(io, nastring)
            end
            if j < p
		print(io, separator)
            else
                print(io, '\n')
            end
        end
    end
    return
end

function printtable(df::AbstractDataFrame;
                    header::Bool = true,
                    separator::Char = ',',
                    quotemark::Char = '"',
                    nastring::AbstractString = "NA")
    printtable(STDOUT,
               df,
               header = header,
               separator = separator,
               quotemark = quotemark,
               nastring = nastring)
    return
end

"""
Write data to a tabular-file format (CSV, TSV, ...)

```julia
writetable(filename, df, [keyword options])
```

### Arguments

* `filename::AbstractString` : the filename to be created
* `df::AbstractDataFrame` : the AbstractDataFrame to be written

### Keyword Arguments

* `separator::Char` -- The separator character that you would like to use. Defaults to the output of `getseparator(filename)`, which uses commas for files that end in `.csv`, tabs for files that end in `.tsv` and a single space for files that end in `.wsv`.
* `quotemark::Char` -- The character used to delimit string fields. Defaults to `'"'`.
* `header::Bool` -- Should the file contain a header that specifies the column names from `df`. Defaults to `true`.
* `nastring::AbstractString` -- What to write in place of missing data. Defaults to `"NA"`.

### Result

* `::DataFrame`

### Examples

```julia
df = DataFrame(A = 1:10)
writetable("output.csv", df)
writetable("output.dat", df, separator = ',', header = false)
writetable("output.dat", df, quotemark = '\'', separator = ',')
writetable("output.dat", df, header = false)
```
"""
function writetable(filename::AbstractString,
                    df::AbstractDataFrame;
                    header::Bool = true,
                    separator::Char = getseparator(filename),
                    quotemark::Char = '"',
                    nastring::AbstractString = "NA",
                    append::Bool = false)

    if endswith(filename, ".bz") || endswith(filename, ".bz2")
        throw(ArgumentError("BZip2 compression not yet implemented"))
    end

    if append && isfile(filename) && filesize(filename) > 0
        file_df = readtable(filename, header = false, nrows = 1)

        # Check if number of columns matches
        if size(file_df, 2) != size(df, 2)
            throw(DimensionMismatch("Number of columns differ between file and DataFrame"))
        end

        # When 'append'-ing to a nonempty file,
        # 'header' triggers a check for matching colnames
        if header
            if any(i -> symbol(file_df[1, i]) != index(df)[i], 1:size(df, 2))
                throw(KeyError("Column names don't match names in file"))
            end

            header = false
        end
    end

    openfunc = endswith(filename, ".gz") ? gzopen : open

    openfunc(filename, append ? "a" : "w") do io
        printtable(io,
                   df,
                   header = header,
                   separator = separator,
                   quotemark = quotemark,
                   nastring = nastring)
    end

    return
end

##############################################################################
#
# HTML output
#
##############################################################################

function html_escape(cell::AbstractString)
    cell = replace(cell, "&", "&")
    cell = replace(cell, "<", "<")
    cell = replace(cell, ">", ">")
    return cell
end

function Base.writemime(io::IO,
                        ::MIME"text/html",
                        df::AbstractDataFrame)
    n = size(df, 1)
    cnames = _names(df)
    write(io, "
")
    write(io, "
")
    write(io, "	")
    for column_name in cnames
        write(io, "$column_name	")
    end
    write(io, "
")
    tty_rows, tty_cols = Base.tty_size()
    mxrow = min(n,tty_rows)
    for row in 1:mxrow
        write(io, "
")
        write(io, "$row	")
        for column_name in cnames
            cell = string(df[row, column_name])
            write(io, "$(html_escape(cell))	")
        end
        write(io, "
")
    end
    if n > mxrow
        write(io, "
")
        write(io, "⋮	")
        for column_name in cnames
            write(io, "⋮	")
        end
        write(io, "
")
    end
    write(io, "
")
end

##############################################################################
#
# MIME
#
##############################################################################

function Base.writemime(io::IO,
                        ::MIME"text/csv",
                        df::AbstractDataFrame)
    printtable(io, df, true, ',')
end

function Base.writemime(io::IO,
                        ::MIME"text/tab-separated-values",
                        df::AbstractDataFrame)
    printtable(io, df, true, '\t')
end


text:
export kahuna_read, kahuna_write, load_dm3, load_dm4, load_hdf5, load_mat, load_mib, load_toml, load_jld2
export save_hdf5, save_toml, save_jld2
export firstheader, type2dict, mib2h5

function kahuna_read(filepath::AbstractString, args...; kwargs...)
    file_extension = splitext(filepath)[2]
    if file_extension == ".dm3"
        output = load_dm3(filepath, args...; kwargs...)
    elseif file_extension == ".dm4"
        output = load_dm4(filepath, args...; kwargs...)
    elseif file_extension == ".hdf5" || file_extension == ".h5"
        output = load_hdf5(filepath, args...; kwargs...)
    elseif file_extension == ".mat"
        output = load_mat(filepath, args...; kwargs...)
    elseif file_extension == ".mib"
        output = load_mib(filepath, args...; kwargs...)
    elseif file_extension == ".toml"
        output = load_toml(filepath, args...; kwargs...)
    elseif file_extension == ".jld2"
        output = load_jld2(filepath, args...; kwargs...)
    else
        output = FileIO.load(filepath)
    end
    return output
end

function kahuna_write(filepath::AbstractString, parent::AbstractString, name::AbstractString, data, args...; kwargs...)
    file_extension = splitext(filepath)[2]
    if file_extension == ".hdf5" || file_extension == ".h5"
        save_hdf5(filepath, parent, name, data, args...; kwargs...)
    elseif file_extension == ".toml"
        save_toml(filepath, args...; kwargs...)
    elseif file_extension == ".jld2"
        save_jld2(filepath, args...; kwargs...)
    else
        FileIO.save(filepath, data)
    end
end

function kahuna_write(filepath::AbstractString, mib_images::Vector{Array{T, 2}} where T <: Union{UInt8, UInt16, UInt32, UInt64}, mib_headers::Vector{MIBHeader}, args...; kwargs...)
    file_extension = splitext(filepath)[2]
    if file_extension == ".hdf5" || file_extension == ".h5"
        save_hdf5(filepath::AbstractString, mib_images::Vector{Array{T, 2}} where T <: Union{UInt8, UInt16, UInt32, UInt64}, mib_headers::Vector{MIBHeader}, args...; kwargs...)
    elseif file_extension == ".jld2"
        @save filepath mib_images mib_headers
    else
        disp("Data not saved.")
    end
end

function load_dm3(filepath)

end

function load_dm4(filepath)

end

function load_hdf5(filepath)

end

function load_mat(filepath::AbstractString, varname::AbstractString)
    file = matopen(filepath)
    if exists(file, varname)
        output = read(file, varname)
    end
    close(file)

    return output
end

function load_mat(filepath::AbstractString; mode="all", kwargs...)
    if mode == "all"
        output = matread(filepath)
    elseif mode == "list"
        file = matopen(filepath)
        output = names(file)
        close(file)
    end
    return output
end

function load_mib(filepath::AbstractString; kwargs...)
    first_header = firstheader(filepath)
    images, headers = read_mib(filepath, first_header; kwargs...)
end

function load_toml(filepath)

end

function load_jld2(filepath)

end

function save_hdf5(filepath::AbstractString, parent::AbstractString, name::AbstractString, data, args...; kwargs...)
    fid = h5open(filepath, "w"; kwargs...)
        # d_write(parent, name, data)
    close(fid)
end

# function save_hdf5(filepath::AbstractString, mib_images::Vector{Array{Integer, 2}}, mib_headers::Vector{AbstractMIBHeader}, args...; range=[1,typemax(Int)], kwargs...)
function save_hdf5(filepath::AbstractString, mib_images::Vector{Array{T, 2}} where T <: Union{UInt8, UInt16, UInt32, UInt64}, mib_headers::Vector{MIBHeader}; range=[1,typemax(Int)], kwargs...)
    for i = 1:length(mib_images)
        if range[1] <= mib_headers[i].id <= range[2]
            h5write(filepath, string("image_" , lpad(i, 8, "0")), mib_images[i])
            h5writeattr(filepath, string("image_" , lpad(i, 8, "0")), type2dict(mib_headers[i]))
        end
    end
end

function save_toml(filepath)

end

function save_jld2(filepath)

end

function read_mib(filepath::AbstractString, first_header::AbstractMIBHeader; range=[1,typemax(Int)])
    offset = first_header.offset
    type = first_header.data_type
    dims = first_header.dims
    raw = first_header.raw
    image_bit_depth = first_header.image_bit_depth

    fid = open(filepath, "r")
    headers = Vector{MIBHeader}()
    if raw
        depth_dict = Dict(1 => UInt8, 6 => UInt8, 12 => UInt16,
                          24 => UInt32, 48 => UInt64)
        type = depth_dict[image_bit_depth]
    end
        buffer = Array{type}(undef, dims[1], dims[2])
        images = Vector{Array{type, 2}}()

    n = 0
    while eof(fid) == false && n < range[2]
            header_string = read(fid, offset)
            read!(fid, buffer)
            n += 1
            if n >= range[1]
                push!(headers, make_mibheader(String(header_string); id=n))
                push!(images, hton.(buffer))
            end
    end
    close(fid)
    return images, headers
end

function firstheader(filepath)
    fid = open(filepath)
    trial = split(String(read(fid, 384)), ",")
    offset = parse(Int, trial[3])
    seekstart(fid)
    header_string = String(read(fid, offset))
    close(fid)
    first_header = make_mibheader(header_string; id=1)
    return first_header
end

function make_mibheader(header_string::AbstractString; id=0)
    header = split(header_string, ",")
    offset = parse(Int, header[3])
    nchip = parse(Int, header[4])
    dims = parse.(Int, header[5:6])
    type_dict = Dict("U1" => UInt8, "U8" => UInt8, "U08" => UInt8, "U16" => UInt16,
                     "U32" => UInt32, "U64" => UInt64, "R64" => UInt64)
    data_type = type_dict[header[7]]
    chip_dims = parse.(Int, split(lstrip(header[8]), "x"))
    time = DateTime(header[10][1:end-3], "y-m-d H:M:S.s")
    exposure_s = parse(Float64, header[11])
    image_bit_depth = parse(Int, header[end-1])
    raw = header[7] == "R64"
    MIBHeader(id,
              offset,
              nchip,
              dims,
              data_type,
              chip_dims,
              time,
              exposure_s,
              image_bit_depth,
              raw)
end

function type2dict(mibtype::AbstractMIBHeader)
    Dict(string(fn) => string(getfield(mibtype, fn)) for fn ∈ fieldnames(typeof(mibtype)))
end

function mib2h5(load_filepath::AbstractString, save_filepath::AbstractString; range=[1,typemax(Int)])
    mib_images, mib_headers = kahuna_read(load_filepath; range=[1,typemax(Int)])
    kahuna_write(save_filepath, mib_images, mib_headers)
end

function mib2h5(load_filepath::AbstractString)
    save_filepath = string(splitext(load_filepath)[1], ".h5")
    mib2h5(load_filepath, save_filepath)
end


text:
# type definitions

export AbstractQC, AbstractQuState, AbstractQuCircuit

# abstracts
abstract AbstractQC{N}
abstract AbstractQuState{T,N}
abstract AbstractQuCircuit{N} <: AbstractQC{N}

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
@ms include("solvers/mip_discret.jl")

function main_dmip(args)
    println("="^70)
    println("Début de l'action dmip")
    inst = Instance(args[:infile])
    sv = MipDiscretSolver(inst)
    solve(sv)

    bestsol = sv.bestsol
    print_sol(bestsol)
    print("Création du fichier \"$(guess_solname(bestsol))\"... ")
    write(bestsol)
    println("FAIT !")

    println("Fin de l'action dmip")
end

main_dmip(Args.args)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/ModelReduction.jl/blob/master/LICENSE

using ModelReduction
using Base.Test

datadir = joinpath(Pkg.dir("ModelReduction"), "test",
                   first(splitext(basename(@__FILE__))))

# Model was constructed using the following code, which is now commented
# to break dependency to JuliaFEM. Assembled global matrices K and M are
# saved to binary format after first test, and can be loaded using JLD.
#
# using JuliaFEM
# using JuliaFEM.Preprocess
# using JuliaFEM.Postprocess
# using JuliaFEM.Abaqus: create_surface_elements
# using Logging
# Logging.configure(level=DEBUG)
#
# # read mesh
# mesh = abaqus_read_mesh(joinpath(datadir, "testmodel.inp"))
# info("element sets = ", collect(keys(mesh.element_sets)))
# info("surface sets = ", collect(keys(mesh.surface_sets)))
#
# # create two field problems with different material properties
# cylinder = Problem(Elasticity, "CYLINDER", 3)
# cylinder.elements = create_elements(mesh, "CYLINDER")
# update!(cylinder.elements, "youngs modulus", 165.0)
# update!(cylinder.elements, "poissons ratio", 0.275)
# update!(cylinder.elements, "density", 7.10E-9)
# block = Problem(Elasticity, "BLOCK", 3)
# block.elements = create_elements(mesh, "BLOCK")
# update!(block.elements, "youngs modulus", 208.0)
# update!(block.elements, "poissons ratio", 0.30)
# update!(block.elements, "density", 7.80E-9)
#
# # create boundary condition from node set
# fixed = Problem(Dirichlet, "fixed", 3, "displacement")
# fixed_nodes = mesh.node_sets[:FIXED]
# fixed.elements = [Element(Poi1, [nid]) for nid in fixed_nodes]
# update!(fixed.elements, "geometry", mesh.nodes)
# update!(fixed.elements, "displacement 1", 0.0)
# update!(fixed.elements, "displacement 2", 0.0)
# update!(fixed.elements, "displacement 3", 0.0)
#
# """ A helper function to create tie contact. """
# function create_interface(mesh::Mesh, slave::String, master::String)
#     interface = Problem(Mortar, "tie contact", 3, "displacement")
#     interface.properties.dual_basis = false
#     slave_elements = create_surface_elements(mesh, slave)
#     master_elements = create_surface_elements(mesh, master)
#     update!(slave_elements, "master elements", master_elements)
#     interface.elements = [slave_elements; master_elements]
#     return interface
# end
#
# # call helper function to create tie contacts
# tie = create_interface(mesh, "CYLINDER_TO_BLOCK", "BLOCK_TO_CYLINDER")
#
# # solve eigenvalues:
# analysis = Analysis(Modal)
# add_problems!(analysis, [cylinder, block, fixed, tie])
#

# skip assembly for now ...
using JLD
data = load(joinpath(datadir, "data.jld"))

cb = Analysis(CraigBampton)
# cb.properties.r_nodes = r_nodes = collect(mesh.node_sets[:BORDER])
# cb.properties.l_nodes = setdiff(keys(mesh.nodes), r_nodes)
# add_problems!(cb, [cylinder, block, fixed, tie])
cb.properties.K = data["K"]
cb.properties.M = data["M"]
cb.properties.r_nodes = data["r_nodes"]
cb.properties.l_nodes = data["l_nodes"]

run!(cb)

info("Calculate 10 smallest eigenvalues of the reduced system")
t0 = time()
w, X = eigs(cb.properties.K_red, cb.properties.M_red; which=:SM, nev=10)
dt = round(time() - t0, 2)
info("Calculated eigenvalues in $dt seconds")

freq_reduced = sqrt.(real(w))/(2*pi)
freq_original = [1937.04, 1938.56, 3612.25, 4240.52, 4254.46, 6036.87, 6224.13, 8835.43, 8876.81, 11716.6]
info("Eigenvalues of the original system [Hz]: ", freq_original)
info("Eigenvalues of the reduced system [Hz]: ", freq_reduced)

@test isapprox(freq_original, freq_reduced; rtol=0.051)