text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.4
#     language: julia
#     name: julia-1.6
# ---

# %%
# https://twitter.com/iitenki_moruten/status/1434070751884177415

using Plots
using StatsBase
using Random
using Printf
using Distributions
using BenchmarkTools
using QuadGK
using SpecialFunctions
using Zygote

# Generating Gaussian distribution (1-dim) by using HMC method
# 改良後

v(x) = 0.5*x^2
u(x) = x^2*(x-3)^2/2

function S(x,n_size)
    sum(u(x[i]) for i=1:n_size)
end

function dSdx(x,k)
    return u'(x)
end

function H(x,p,n_size)
    K = sum(p[i]^2 for i=1:n_size)
    return S(x,n_size) + 0.5*K
end

#１変数更新のHMC法
function HMC_1dim(x,parameter)
        n_size,n_tau,τ = parameter
        #初期化
        p_ini = rand(Normal(0,1),n_size)
        p_fin = zeros(n_size)
        x_ini = zeros(n_size)
        x_fin = zeros(n_size)

        #初期配位を保存
        for i=1:n_size x_ini[i] = x[i] end

        #初期配位のハミルトニアンを計算
        H_ini = H(x_ini,p_ini,n_size)

        #リープフロッグで時間発展
        leapfrog!(x,p_ini,parameter)

        #時間発展後の配位を保存
        for i=1:n_size x_fin[i] = x[i]; p_fin[i] = p_ini[i] end

        #時間発展後のハミルトニアンを計算
        H_fin = H(x_fin,p_fin,n_size)

        #Metropolis check
        r=rand()
        if r<exp(H_ini - H_fin)
            for i=1:n_size x[i] = x_fin[i] end
        else
            for i=1:n_size x[i] = x_ini[i] end
        end
end

#x,pをリープフロッグ法で更新する．
function leapfrog!(x,p,parameter)
    n_size,n_tau,τ = parameter
    Δτ  = τ/n_tau
    p0  = 0.0
    p1  = 0.0
    x05 = 0.0
    x15 = 0.0
    for k=1:n_size
        #initial setting : τ0
        p0  = p[k]
        x05 = x[k] + p0 * 0.5*Δτ

        #(n-1)step : τ0 -> τn-1
        for i=1:n_tau-1
            p1  = p0 - dSdx(x05,k)*Δτ
            x15 = x05 + p1*Δτ
            p0  = p1
            x05 = x15
        end

        #final step : τn-1 -> τn
        p_end = p1 - dSdx(x15,k)*Δτ
        x_end = x15 + p_end*0.5*Δτ

        x[k]  = x_end
        p[k]  = p_end
    end
end

function gen_configs(Ω,n_con,parameter)
    n_size,n_tau,τ = parameter
    x = zeros(n_size)
    for j=1:n_con
        HMC_1dim(x,parameter)
        Ω[j] = x[1]
    end
end

##########################################################

"""q = [n_size,n_tau,τ]"""

#パラメータベクトル
n_config = 1_000_000          # # of configurations
n_size   = 1
n_tau    = 100
τ        = 5
q 　　　  = [n_size,n_tau,τ]
Ω        = zeros(n_config)

#result = @benchmark gen_configs(n_config,q)
@time gen_configs(Ω,n_config,q)

##########################################################

###plot###
gauss(x) = exp(-0.5*x^2)/sqrt(2*pi)
pdf_u(x) = exp(-u(x))
Z0,error = quadgk(pdf_u,-Inf,Inf)
normed_u(x) = pdf_u(x)/Z0

xr = range(-3,6,length=1000)

plot(Ω,st=:histogram,nbins=50,norm=:pdf,alpha=0.3,label="HMC's sample",legend=:topleft)
plot!(xr,normed_u.(xr),label="pdf_ana")

# %%
##########################################################

"""q = [n_size,n_tau,τ]"""

#パラメータベクトル
n_config = 1_000_000          # # of configurations
n_size   = 1
n_tau    = 100
τ        = 5
q 　　　  = [n_size,n_tau,τ]
Ω        = zeros(n_config)

#result = @benchmark gen_configs(n_config,q)
@time gen_configs(Ω,n_config,q)

##########################################################

###plot###
gauss(x) = exp(-0.5*x^2)/sqrt(2*pi)
pdf_u(x) = exp(-u(x))
Z0,error = quadgk(pdf_u,-Inf,Inf)
normed_u(x) = pdf_u(x)/Z0

xr = range(-3,6,length=1000)

plot(Ω,st=:histogram,nbins=50,norm=:pdf,alpha=0.3,label="HMC's sample",legend=:topleft)
plot!(xr,normed_u.(xr),label="pdf_ana")

# %% [markdown]
# 以下のコードでの主な変更点
#
# * `r < exp(H_ini - H_fin)` → `r ≤ min(1, exp(H_ini - H_fin))` (Metropolis-Hastings化)
# * `n_size` を函数内から削除 (必要なら `n_size = length(x)` とできる)
# * `parameter` に作業用配列 `p_ini`, `p_fin`, `x_ini`, `x_fin` を含めて、メモリアロケーションを削減。
# * `parameter` にポテンシャル函数 `u` も含め、函数 `H` などに渡すようにした。
# * Zygote.jl ではなく、ForwardDiff.jl を使うようにして、初回実行時の遅延の短縮。
#
# References
#
# * https://docs.julialang.org/en/v1/manual/performance-tips/#Pre-allocating-outputs
# * https://docs.julialang.org/en/v1/manual/performance-tips/#More-dots:-Fuse-vectorized-operations
# * https://github.com/genkuroki/public/blob/main/0018/HMC%20leapfrog.ipynb

# %%
using Plots
# using StatsBase
# using Random
# using Printf
using Distributions
using BenchmarkTools
using QuadGK
# using SpecialFunctions

module Revised

using Random
using Distributions
using ForwardDiff

# Generating Gaussian distribution (1-dim) by using HMC method
# 改良後

S(x, parameter) = ((u,) = parameter; sum(u, x))
dSdx(x, k, parameter) = ((u,) = parameter; ForwardDiff.derivative(u, x))
square(x) = x^2
H(x, p, parameter) = S(x, parameter) + 0.5sum(square, p)

#１変数更新のHMC法
function HMC_1dim!(x, parameter)
        _, n_tau, τ, p_ini, p_fin, x_ini, x_fin = parameter

        #初期化
        rand!(Normal(0,1), p_ini)

        #初期配位を保存
        x_ini .= x

        #初期配位のハミルトニアンを計算
        H_ini = H(x_ini, p_ini, parameter)

        #リープフロッグで時間発展
        leapfrog!(x, p_ini, parameter)

        #時間発展後の配位を保存
        x_fin .= x
        p_fin .= p_ini

        #時間発展後のハミルトニアンを計算
        H_fin = H(x_fin, p_fin, parameter)

        #Metropolis check
        r = rand()
        if r ≤ min(1, exp(H_ini - H_fin))
            x .= x_fin
        else
            x .= x_ini
        end
end

#x,pをリープフロッグ法で更新する．
function leapfrog!(x, p, parameter)
    n_size = length(x)
    _, n_tau, τ = parameter
    Δτ  = τ/n_tau
    for k=1:n_size
        #initial setting : τ0
        p0  = p[k]
        x05 = x[k] + p0 * 0.5*Δτ

        #(n-1)step : τ0 -> τn-1
        local p1, x15
        for i=1:n_tau-1
            p1  = p0 - dSdx(x05, k, parameter)*Δτ
            x15 = x05 + p1*Δτ
            p0  = p1
            x05 = x15
        end

        #final step : τn-1 -> τn
        p_end = p1 - dSdx(x15, k, parameter)*Δτ
        x_end = x15 + p_end*0.5*Δτ

        x[k]  = x_end
        p[k]  = p_end
    end
end

function gen_configs(Ω, x, parameter)
    n_tau, τ = parameter
    for j in eachindex(Ω)
        HMC_1dim!(x, parameter)
        Ω[j] = x[1]
    end
end

end

v(x) = 0.5x^2
u(x) = x^2*(x-3)^2/2

##########################################################

n_config = 10^6          # # of configurations
n_size   = 1
n_tau    = 100
τ        = 5

Ω        = zeros(n_config)
x        = zeros(n_size)
p_ini    = similar(x)
p_fin    = similar(x)
x_ini    = similar(x)
x_fin    = similar(x)

"""q = (u, n_tau, τ, p_ini, p_fin, x_ini, x_fin)"""
q 　　　  = (u, n_tau, τ, p_ini, p_fin, x_ini, x_fin)

#result = @benchmark gen_configs(n_config,q)
@time Revised.gen_configs(Ω, x, q)

##########################################################

###plot###
gauss(x) = exp(-0.5*x^2)/sqrt(2*pi)
pdf_u(x) = exp(-u(x))
Z0,error = quadgk(pdf_u,-Inf,Inf)
normed_u(x) = pdf_u(x)/Z0

xr = range(-3,6,length=1000)

plot(Ω,st=:histogram,nbins=50,norm=:pdf,alpha=0.3,label="HMC's sample",legend=:topleft)
plot!(xr,normed_u.(xr),label="pdf_ana")

# %%
##########################################################

n_config = 10^6          # # of configurations
n_size   = 1
n_tau    = 100
τ        = 5

Ω        = zeros(n_config)
x        = zeros(n_size)
p_ini    = similar(x)
p_fin    = similar(x)
x_ini    = similar(x)
x_fin    = similar(x)

"""q = (u, n_tau, τ, p_ini, p_fin, x_ini, x_fin)"""
q 　　　  = (u, n_tau, τ, p_ini, p_fin, x_ini, x_fin)

#result = @benchmark gen_configs(n_config,q)
@time Revised.gen_configs(Ω, x, q)

##########################################################

###plot###
gauss(x) = exp(-0.5*x^2)/sqrt(2*pi)
pdf_u(x) = exp(-u(x))
Z0,error = quadgk(pdf_u,-Inf,Inf)
normed_u(x) = pdf_u(x)/Z0

xr = range(-3,6,length=1000)

plot(Ω,st=:histogram,nbins=50,norm=:pdf,alpha=0.3,label="HMC's sample",legend=:topleft)
plot!(xr,normed_u.(xr),label="pdf_ana")

# %%

max_stars_repo_path: 0019/HMC leapfrog.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0019/HMC leapfrog.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0019/HMC leapfrog.jl

max_forks_repo_name: genkuroki/public


text:
function print_synthetic_fasta(path::String,line_length::Integer=80)
	header=">CM002885.2.1 BioBackgroundModels Synthetic chromosome for tests\n"
	write(path, header, format_lines(generate_synthetic_seq(),line_length))
end

function print_synthetic_index(path::String)
	write(path, "CM002885.2.1	1000	5	80	81\n")
end

function print_synthetic_position(path::String)
	write(path, ">1 start 501 end 641 smt_pos 111 smt_value 205.11821 fuzziness_score 43.53698\n",generate_synthetic_seq()[501:641],"\n")
end

function print_synthetic_gff(path::String)
gff_txt = "##gff-version 3
##sequence-region   1 1 1000
#!genome-build BGHMM Synthetic
#!genome-version Synth1.0
#!genome-date 2019-08
#!genome-build-accession NCBI:FAKE
#!genebuild-last-updated 2019-08
1	Ensembl	chromosome	1	1000	.	.	.	ID=chromosome:1;Alias=CM002885.2,NC_000001.1
###
1	ensembl_havana	gene	501	1000	.	-	.	ID=gene:ENSDARG00000000001;Name=fake01;biotype=protein_coding;description=fake gene [Source:FAKESRC];gene_id=ENSDARG00000000001;logic_name=ensembl_havana_gene;version=1
1	ensembl_havana	mRNA	501	1000	.	-	.	ID=transcript:ENSDART00000000001;Parent=gene:ENSDARG00000000001;Name=fake01-203;biotype=protein_coding;transcript_id=ENSDART00000000001;version=1
1	ensembl_havana	three_prime_UTR	501	510	.	-	.	Parent=transcript:ENSDART00000000001
1	ensembl_havana	exon	501	570	.	+	.	Parent=transcript:ENSDART00000000001;Name=ENSDARE00000000001;constitutive=0;ensembl_end_phase=0;ensembl_phase=0;exon_id=ENSDARE00000000001;rank=5;version=1
1	ensembl_havana	CDS	511	570	.	+	0	ID=CDS:ENSDARP00000000001;Parent=transcript:ENSDART00000000001;protein_id=ENSDARP00000000001
1	ensembl_havana	exon	601	660	.	-	.	Parent=transcript:ENSDART00000000001;Name=ENSDARE00000000002;constitutive=0;ensembl_end_phase=0;ensembl_phase=0;exon_id=ENSDARE00000000002;rank=4;version=1
1	ensembl_havana	CDS	601	660	.	-	1	ID=CDS:ENSDARP00000000001;Parent=transcript:ENSDART00000000001;protein_id=ENSDARP00000000001
1	ensembl_havana	exon	701	760	.	-	.	Parent=transcript:ENSDART00000000001;Name=ENSDARE00000000003;constitutive=0;ensembl_end_phase=0;ensembl_phase=0;exon_id=ENSDARE00000000003;rank=3;version=1
1	ensembl_havana	CDS	701	760	.	-	0	ID=CDS:ENSDARP00000000001;Parent=transcript:ENSDART00000000001;protein_id=ENSDARP00000000001
1	ensembl_havana	exon	801	860	.	-	.	Parent=transcript:ENSDART00000000001;Name=ENSDARE00000000004;constitutive=0;ensembl_end_phase=0;ensembl_phase=0;exon_id=ENSDARE00000000004;rank=2;version=1
1	ensembl_havana	CDS	801	860	.	-	0	ID=CDS:ENSDARP00000000001;Parent=transcript:ENSDART00000000001;protein_id=ENSDARP00000000001
1	ensembl_havana	exon	901	975	.	-	.	Parent=transcript:ENSDART00000000001;Name=ENSDARE00000000005;constitutive=0;ensembl_end_phase=0;ensembl_phase=0;exon_id=ENSDARE00000000005;rank=1;version=1
1	ensembl_havana	CDS	901	960	.	-	1	ID=CDS:ENSDARP00000000001;Parent=transcript:ENSDART00000000001;protein_id=ENSDARP00000000001
1	ensembl_havana	five_prime_UTR	961	975	.	-	.	Parent=transcript:ENSDART00000000001
###\n"
	write(path,gff_txt)
end

function generate_synthetic_seq(gene_start::Integer=501, UTR3L::Integer=10, exon_length::Integer=60, intron_length::Integer=40, no_exons::Integer=5, UTR5L::Integer=25, OAL::Integer=1000, (iseq,pseq,eseq)::Tuple{String,String,String}=("AT","CG","CAT"); verbose::Bool=false)
	@assert length(iseq) == 2
	@assert length(pseq) == 2
	@assert length(eseq) == 3
	@assert mod(gene_start-1,2) == 0
	@assert mod(exon_length,3) == 0
	@assert mod(intron_length,2) == 0
	gene_length=gene_start+UTR3L+(no_exons*exon_length)+((no_exons-1)*intron_length)+UTR5L
	@assert gene_length <=OAL
	seq=""
	for i in 1:((gene_start-1) / 2)
		seq *= iseq
	end
	verbose && @info "Intergenic length $(length(seq))"
	for i in 1:UTR3L / 2
		seq *= pseq
	end
	verbose && @info "Added 3'UTR... $(length(seq))"
	for ex in 1:no_exons-1
		for i in 1:(exon_length/3)
			seq*=eseq
		end
		ex == 1 ? ilength = 30 : ilength = 40
		for i in 1:(ilength/2)
			seq*=pseq
		end
	end
	for i in 1:(exon_length/3)
		seq*=eseq
	end
	verbose && @info "Added exons and introns... $(length(seq))"
	for i in 1:UTR5L / 2
		seq *= pseq
	end
	verbose && @info "Added 5'UTR... $(length(seq))"

	for i in 1:((OAL-length(seq))/2)
		seq*=iseq
	end

	verbose && @info "Generated synthetic genome sequence of length $(length(seq))"

	return seq
end

  function format_lines(seq::String, line_length::Integer)
	a=join((SubString(seq,i,min(i+line_length-1,length(seq))) for i=1:line_length:length(seq)),'\n')
	return(a)
  end

ext: jl


text:
#!/usr/bin/env julia
# > julia -O3 geoDistMP.jl

#---------------------------------------------------------------------
# [北緯]	[東経]	[その他]
#---------------------------------------------------------------------
# 十進法 ddd.d...
#---------------------------------------------------------------------
const Data1 = """
// (サンプル)
//	35.685187	139.752274	皇居
//	24.449582	122.934340	日本最西端之地

// (参考)
//	https://www.330k.info/essay/precision-of-mathematica-geodistance-and-geographical-distance-formula/
//	https://www.330k.info/essay/precision-of-mathematica-geodistance-and-geographical-distance-formula/makurazaki_wakkanai_reduced.gpx
31.27216	130.29956	枕崎駅
31.33976	130.27063	↓
31.38294	130.31228	↓
31.45631	130.34415	↓
31.58710	130.34180	↓
31.66086	130.32140	↓
31.71886	130.27254	↓
31.75990	130.30652	↓
31.81231	130.30765	↓
31.84331	130.27986	↓
31.85158	130.23748	↓
31.86527	130.22490	↓
31.92860	130.22425	↓
32.02740	130.19569	↓
32.05986	130.22292	↓
32.11650	130.33952	↓
32.17422	130.38257	↓
32.19819	130.38227	↓
32.25715	130.50983	↓
32.33472	130.48537	↓
32.36339	130.50789	↓
32.37484	130.53627	↓
32.42086	130.55128	↓
32.45672	130.59943	↓
32.50328	130.61136	↓
32.50747	130.64885	↓
32.58873	130.70275	↓
32.65366	130.69925	↓
32.68998	130.67069	↓
32.82402	130.72266	↓
32.89552	130.68972	↓
32.96968	130.71602	↓
33.00075	130.68708	↓
33.02791	130.68218	↓
33.10306	130.70025	↓
33.11602	130.68786	↓
33.11307	130.65209	↓
33.13460	130.62086	↓
33.19131	130.61054	↓
33.21966	130.55923	↓
33.31397	130.52229	↓
33.44212	130.53532	↓
33.45249	130.56446	↓
33.50689	130.60107	↓
33.51256	130.63301	↓
33.55673	130.64264	↓
33.59143	130.67364	↓
33.64598	130.67801	↓
33.66647	130.71552	↓
33.75130	130.75641	↓
33.79189	130.74460	↓
33.86573	130.76159	↓
33.86768	130.81787	↓
33.88142	130.83614	↓
33.87258	130.88045	↓
33.89303	130.92202	↓
33.94799	130.96864	↓
33.97642	130.94679	↓
33.99221	130.95232	↓
34.00376	130.99296	↓
34.06037	131.03391	↓
34.04489	131.06879	↓
34.05162	131.19582	↓
34.03629	131.25252	↓
34.04143	131.27890	↓
34.06003	131.29025	↓
34.05707	131.36656	↓
34.08422	131.39656	↓
34.10589	131.39556	↓
34.15606	131.43314	↓
34.20717	131.51493	↓
34.25024	131.52257	↓
34.28446	131.55714	↓
34.34519	131.63228	↓
34.35404	131.67313	↓
34.41413	131.74881	↓
34.44155	131.77275	↓
34.47418	131.77728	↓
34.48511	131.80555	↓
34.54401	131.83453	↓
34.55209	131.83068	↓
34.54495	131.81727	↓
34.56998	131.79227	↓
34.62423	131.80686	↓
34.63737	131.79758	↓
34.64592	131.81214	↓
34.67804	131.81309	↓
34.70833	131.86113	↓
34.74774	131.87338	↓
34.77557	131.96065	↓
34.82374	131.98066	↓
34.86535	132.02724	↓
34.89566	132.10535	↓
34.93643	132.11146	↓
34.99852	132.19486	↓
35.01339	132.23355	↓
35.03471	132.24472	↓
35.03723	132.28156	↓
35.06997	132.32326	↓
35.07125	132.34562	↓
35.18688	132.42729	↓
35.20403	132.49928	↓
35.23630	132.53112	↓
35.24897	132.56902	↓
35.27473	132.58686	↓
35.29587	132.64745	↓
35.33366	132.68066	↓
35.34675	132.72335	↓
35.36931	132.73871	↓
35.37298	132.79000	↓
35.39050	132.81830	↓
35.45009	133.08633	↓
35.42799	133.13184	↓
35.43300	133.24658	↓
35.42482	133.28367	↓
35.40838	133.29836	↓
35.40927	133.32596	↓
35.42741	133.41244	↓
35.46674	133.44513	↓
35.50942	133.52286	↓
35.50877	133.63961	↓
35.48485	133.69113	↓
35.50181	133.72706	↓
35.49587	133.89496	↓
35.51766	133.96873	↓
35.51196	134.00909	↓
35.52838	134.14519	↓
35.52144	134.23435	↓
35.56746	134.29629	↓
35.55688	134.36694	↓
35.52550	134.39883	↓
35.52472	134.42653	↓
35.56568	134.48173	↓
35.49301	134.54081	↓
35.49372	134.56549	↓
35.47117	134.59804	↓
35.41192	134.60040	↓
35.37498	134.64070	↓
35.39786	134.76937	↓
35.37299	134.78917	↓
35.31739	134.88940	↓
35.33466	134.91473	↓
35.31991	134.98540	↓
35.33376	135.09495	↓
35.39405	135.15019	↓
35.39259	135.18012	↓
35.42419	135.20853	↓
35.42966	135.24969	↓
35.45320	135.26569	↓
35.45204	135.34753	↓
35.47645	135.36740	↓
35.48973	135.52341	↓
35.47780	135.65104	↓
35.49061	135.74653	↓
35.46257	135.86101	↓
35.47916	135.87757	↓
35.47875	135.90096	↓
35.58592	135.91037	↓
35.59768	135.92192	↓
35.61310	135.97854	↓
35.63258	136.00405	↓
35.63885	136.08300	↓
35.70332	136.07857	↓
35.72268	136.09934	↓
35.75652	136.10639	↓
35.79934	136.09123	↓
35.84969	136.12614	↓
35.88957	136.19207	↓
35.96105	136.18821	↓
36.01722	136.21089	↓
36.04990	136.24152	↓
36.28876	136.29599	↓
36.34476	136.44836	↓
36.38915	136.49315	↓
36.42353	136.48314	↓
36.47515	136.50073	↓
36.51402	136.57280	↓
36.53943	136.60038	↓
36.56844	136.60161	↓
36.63330	136.70551	↓
36.68492	136.72455	↓
36.67215	136.80055	↓
36.68909	136.86221	↓
36.68596	136.89412	↓
36.76182	137.00545	↓
36.75750	137.06347	↓
36.72641	137.13097	↓
36.73463	137.20904	↓
36.72098	137.25304	↓
36.71967	137.31469	↓
36.78871	137.40655	↓
36.83859	137.42126	↓
36.87137	137.45589	↓
36.91488	137.47441	↓
36.93896	137.51278	↓
36.94798	137.56442	↓
36.97093	137.58841	↓
36.99355	137.70345	↓
37.05869	137.91101	↓
37.12734	138.04818	↓
37.17272	138.09649	↓
37.16073	138.15155	↓
37.16793	138.26336	↓
37.31293	138.42690	↓
37.34953	138.49185	↓
37.36877	138.57680	↓
37.44339	138.63320	↓
37.46300	138.67949	↓
37.56251	138.73300	↓
37.59259	138.77035	↓
37.63742	138.86993	↓
37.77192	138.90211	↓
37.78816	138.92311	↓
37.84504	138.93742	↓
37.93167	139.13773	↓
37.96558	139.24971	↓
37.99055	139.24460	↓
38.07437	139.34672	↓
38.14536	139.40519	↓
38.14568	139.42419	↓
38.18620	139.44732	↓
38.18575	139.46358	↓
38.23906	139.50305	↓
38.26969	139.50492	↓
38.32283	139.55531	↓
38.40567	139.55791	↓
38.43739	139.59159	↓
38.46184	139.57374	↓
38.48853	139.52186	↓
38.50020	139.52157	↓
38.64672	139.59635	↓
38.68556	139.62721	↓
38.71575	139.67726	↓
38.72160	139.71889	↓
38.71171	139.74665	↓
38.76013	139.83237	↓
38.81958	139.83292	↓
38.84014	139.84887	↓
38.88278	139.83798	↓
38.89732	139.85534	↓
38.92997	139.85910	↓
38.95294	139.83711	↓
39.07196	139.88639	↓
39.11331	139.87010	↓
39.17766	139.90869	↓
39.22050	139.90142	↓
39.25122	139.92663	↓
39.27512	139.92511	↓
39.29048	139.97351	↓
39.38573	140.03200	↓
39.42543	140.02743	↓
39.49556	140.04930	↓
39.70086	140.06837	↓
39.72331	140.08769	↓
39.74754	140.06741	↓
39.84763	140.05160	↓
39.88920	140.08955	↓
39.94486	140.11257	↓
39.94495	140.20285	↓
39.95930	140.23087	↓
39.99105	140.24643	↓
40.00333	140.27761	↓
40.03392	140.27518	↓
40.06666	140.29524	↓
40.11973	140.37870	↓
40.17704	140.40089	↓
40.25199	140.38117	↓
40.27930	140.45303	↓
40.27213	140.53108	↓
40.34817	140.60031	↓
40.41502	140.61142	↓
40.42731	140.63640	↓
40.46110	140.63580	↓
40.50639	140.61219	↓
40.52976	140.55497	↓
40.55273	140.55145	↓
40.58602	140.50747	↓
40.64064	140.48979	↓
40.66853	140.51135	↓
40.71051	140.57928	↓
40.76500	140.58829	↓
40.78177	140.62460	↓
40.80128	140.62282	↓
40.83171	140.66296	↓
40.82293	140.74535	↓
40.82983	140.81478	↓
40.84728	140.83660	↓
40.91711	140.87235	↓
40.91844	140.93897	↓
40.93293	140.99256	↓
40.90375	141.06983	↓
40.87306	141.10033	↓
40.87277	141.13462	↓
40.93775	141.20909	↓
41.01692	141.24443	↓
41.09996	141.25151	↓
41.15419	141.28298	↓
41.18053	141.28334	↓
41.27993	141.20937	↓
41.30066	141.22201	↓
41.37187	141.20907	↓
41.39627	141.17397	↓
41.39664	141.15387	↓
41.41286	141.16244	↓
41.46475	141.10580	↓
41.48141	141.06773	↓
41.48537	141.00212	↓
41.52699	140.93774	↓
41.52438	140.89957	↓
41.75257	140.68421	↓
41.78182	140.68459	↓
41.81038	140.71472	↓
41.90217	140.63790	↓
41.96980	140.65151	↓
42.05680	140.59637	↓
42.07894	140.59760	↓
42.20767	140.40268	↓
42.24245	140.29294	↓
42.26427	140.27304	↓
42.30607	140.27232	↓
42.45768	140.33426	↓
42.52100	140.38371	↓
42.57805	140.45371	↓
42.59251	140.45217	↓
42.60211	140.48068	↓
42.59356	140.53436	↓
42.60455	140.56251	↓
42.59114	140.59347	↓
42.60339	140.61638	↓
42.58390	140.72731	↓
42.61234	140.72931	↓
42.61843	140.76694	↓
42.65231	140.79488	↓
42.66851	140.84224	↓
42.69185	140.83422	↓
42.70496	140.85840	↓
42.73792	140.87731	↓
42.75720	140.92583	↓
42.79565	140.93870	↓
42.80135	140.99045	↓
42.84834	141.02581	↓
42.84628	141.08122	↓
42.87149	141.11492	↓
42.87555	141.14656	↓
42.90214	141.11848	↓
42.93042	141.11944	↓
42.96201	141.15935	↓
42.96970	141.20916	↓
42.95558	141.24633	↓
42.96754	141.33122	↓
43.01499	141.34674	↓
43.06054	141.34076	↓
43.07076	141.39598	↓
43.09376	141.42549	↓
43.09500	141.45322	↓
43.11657	141.47082	↓
43.13105	141.50724	↓
43.19031	141.52474	↓
43.21205	141.51812	↓
43.26471	141.59038	↓
43.31566	141.60110	↓
43.36232	141.73306	↓
43.42824	141.81820	↓
43.46746	141.84986	↓
43.54583	141.88487	↓
43.58899	141.86815	↓
43.66067	141.89507	↓
43.73561	141.87605	↓
43.76817	141.88190	↓
43.79619	141.83838	↓
43.84670	141.82205	↓
43.85322	141.78418	↓
43.84608	141.76329	↓
43.88505	141.73457	↓
43.90817	141.69254	↓
43.91968	141.69585	↓
43.91999	141.67583	↓
43.96137	141.64284	↓
44.03005	141.66835	↓
44.30977	141.65418	↓
44.42645	141.74799	↓
44.54555	141.76696	↓
44.55967	141.78596	↓
44.62337	141.80358	↓
44.79147	141.79018	↓
44.89567	141.74657	↓
44.97043	141.76426	↓
45.00096	141.84433	↓
45.05006	141.82385	↓
45.08409	141.78798	↓
45.10466	141.79712	↓
45.20012	141.77631	↓
45.22762	141.74642	↓
45.24091	141.75248	↓
45.28981	141.72318	↓
45.33368	141.72053	↓
45.35275	141.73673	↓
45.39518	141.70627	↓
45.40446	141.67576	↓
45.41726	141.67677	稚内駅
"""
#---------------------------------------------------------------------
# 度分秒 dddmmss.s...
#---------------------------------------------------------------------
const Data2 = """
// (サンプル)
354106.6732	1394508.1864	皇居
242658.4952	1225603.6240	日本最西端之地
"""
#---------------------------------------------------------------------

#-------------------
# 度分秒 => 十進法
#-------------------
function
rtnGeoIBLto10B(
	ddmmss::Float64 # ddmmss.s...
)
	sign = 1

	if ddmmss < 0
		sign = -1
		ddmmss = -ddmmss
	end

	sec = ddmmss % 100
	min = floor(ddmmss / 100) % 100
	deg = floor(ddmmss / 10000)

	return sign * (deg + (min / 60.0) + (sec / 3600.0))::Float64
end

#-------------------------------
# Vincenty法による２点間の距離
#-------------------------------
function
rtnGeoVincentry(
	lat1::Float64, # 開始～緯度
	lng1::Float64, # 開始～経度
	lat2::Float64, # 終了～緯度
	lng2::Float64  # 終了～経度
)
	if lat1 == lat2 && lng1 == lng2
		return (0.0, 0.0)
	end

	## _A = 6378137.0
	_B   = 6356752.314140356    # GRS80
	_F   = 0.003352810681182319 # 1 / 298.257222101
	_RAD = 0.017453292519943295 # π / 180

	latR1 = lat1 * _RAD
	lngR1 = lng1 * _RAD
	latR2 = lat2 * _RAD
	lngR2 = lng2 * _RAD

	f1 = 1 - _F

	omega  = lngR2 - lngR1
	tanU1  = f1 * tan(latR1)
	cosU1  = 1 / sqrt(1 + tanU1 * tanU1)
	sinU1  = tanU1 * cosU1
	tanU2  = f1 * tan(latR2)
	cosU2  = 1 / sqrt(1 + tanU2 * tanU2)
	sinU2  = tanU2 * cosU2

	lamda     = omega
	dLamda    = 0.0
	sinLamda  = 0.0
	cosLamda  = 0.0
	sin2sigma = 0.0
	sinSigma  = 0.0
	cosSigma  = 0.0
	sigma     = 0.0
	sinAlpha  = 0.0
	cos2alpha = 0.0
	cos2sm    = 0.0
	c = 0.0

	iLoop = 0

	while true
		sinLamda = sin(lamda)
		cosLamda = cos(lamda)
		sin2sigma = (cosU2 * sinLamda) * (cosU2 * sinLamda) + (cosU1 * sinU2 - sinU1 * cosU2 * cosLamda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLamda)
		if sin2sigma < 0.0
			return (0.0, 0.0)
		end
		sinSigma = sqrt(sin2sigma)
		cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLamda
		sigma = atan(sinSigma, cosSigma)
		sinAlpha = cosU1 * cosU2 * sinLamda / sinSigma
		cos2alpha = 1 - sinAlpha * sinAlpha
		cos2sm = cosSigma - 2 * sinU1 * sinU2 / cos2alpha
		c = _F / 16 * cos2alpha * (4 + _F * (4 - 3 * cos2alpha))
		dLamda = lamda
		lamda = omega + (1 - c) * _F * sinAlpha * (sigma + c * sinSigma * (cos2sm + c * cosSigma * (-1 + 2 * cos2sm * cos2sm)))

		if abs(lamda - dLamda) <= 1e-12
			break
		end

		iLoop += 1

		# 日本国内であれば５回程度で収束
		if iLoop > 10
			return (-1, -1) # Err
		end
	end

	d2 = cos2alpha * (1 - f1 * f1) / (f1 * f1)
	a = 1 + d2 / 16384 * (4096 + d2 * (-768 + d2 * (320 - 175 * d2)))
	b = d2 / 1024 * (256 + d2 * (-128 + d2 * (74 - 47 * d2)))
	dSigma = b * sinSigma * (cos2sm + b / 4 * (cosSigma * (-1 + 2 * cos2sm * cos2sm) - b / 6 * cos2sm * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2sm * cos2sm)))
	angle = atan(cosU2 * sinLamda, cosU1 * sinU2 - sinU1 * cosU2 * cosLamda) * 57.29577951308232
	dist = _B * a * (sigma - dSigma)

	# 変換
	if angle < 0
		angle += 360.0 # 360度表記
	end
	dist /= 1000.0 # m => km

	return (dist::Float64, angle::Float64)
end

#-------
# Main
#-------
using Printf

const Separater = "\t"

#---------------
# 計算／十進法
#---------------
function
main_Data1()
	iTotalDist = 0.0
	aOld = []

	for _s1 in split(Data1, "\n")
		_s1 = strip(_s1)

		if length(_s1) > 0 && SubString(_s1, 1, 2) != "//"
			as1 = split(_s1, Separater)

			dist, angle = (
				size(aOld)[1] > 0 ?
				rtnGeoVincentry(parse(Float64, aOld[1]), parse(Float64, aOld[2]), parse(Float64, as1[1]), parse(Float64, as1[2])) :
				(0.0, 0.0)
			)

			str = @sprintf("%.6fkm%s%.6f度", dist, Separater, angle)
			for _s2 in as1
				str *= Separater * _s2
			end
			println(str)

			iTotalDist += dist
			aOld = as1
		end
	end

	@printf("%.6fkm\n\n", iTotalDist)
end

#---------------
# 計算／度分秒
#---------------
function
main_Data2()
	iTotalDist = 0.0
	aOld = []

	for _s1 in split(Data2, "\n")
		_s1 = strip(_s1)

		if length(_s1) > 0 && SubString(_s1, 1, 2) != "//"
			as1 = split(_s1, Separater)
			ad1 = []

			for _s2 in as1[1:2]
				push!(ad1, rtnGeoIBLto10B(parse(Float64, _s2)))
			end

			dist, angle = (
				size(aOld)[1] > 0 ?
				rtnGeoVincentry(aOld[1], aOld[2], ad1[1], ad1[2]) :
				(0.0, 0.0)
			)

			str = @sprintf("%.6fkm%s%.6f度", dist, Separater, angle)
			for _s2 in as1
				str *= Separater * _s2
			end
			println(str)

			iTotalDist += dist
			aOld = ad1
		end
	end

	@printf("%.6fkm\n\n", iTotalDist)
end

#-------
# Exec
#-------
main_Data1()
main_Data2()

max_stars_repo_name: iwm-iwama/leaflet-

max_issues_repo_name: iwm-iwama/leaflet-

max_forks_repo_name: iwm-iwama/leaflet-


text:
"Exception type for error thrown by the NetCDF library"
mutable struct NetCDFError <: Exception
    code::Cint
    msg::String
end

"Construct a NetCDFError from the error code"
NetCDFError(code::Cint) = NetCDFError(code, unsafe_string(nc_strerror(code)))

function Base.showerror(io::IO, err::NetCDFError)
    print(io, "NetCDF error code $(err.code):\n\t$(err.msg)")
end

"Check the NetCDF error code, raising an error if nonzero"
function check(code::Cint)
    # zero means success, return
    if code == Cint(0)
        return
    # otherwise throw an error message
    else
        throw(NetCDFError(code))
    end
end

const funext = [
    (Float64, "double", "float64a"),
    (Float32, "float", "float32a"),
    (Int32,   "int", "int32a"),
    (UInt8,   "uchar", "uint8a"),
    (Int8,    "schar", "int8a"),
    (Int16,   "short", "int16a"),
    (Int64,   "longlong", "int64a")]

const ida          = zeros(Int32,1)
const namea        = zeros(UInt8,NC_MAX_NAME+1)
const lengtha      = zeros(Csize_t,1)
const dimida       = zeros(Int32,NC_MAX_VAR_DIMS)
const chunk_sizea  = zeros(Csize_t,NC_MAX_VAR_DIMS)
const storagep     = zeros(Int32,1)
const ndima        = zeros(Int32,1)
const nvara        = zeros(Int32,1)
const ngatta       = zeros(Int32,1)
const nunlimdimida = zeros(Int32,1)
const typea        = zeros(Int32,1)
const natta        = zeros(Int32,1)
const nvals        = zeros(Csize_t,1)
const dima         = zeros(Int32,1)
const varida       = zeros(Int32,1)
const vara         = zeros(Int32,1)
const dumids       = zeros(Int32,NC_MAX_DIMS)
const gstart       = zeros(UInt,NC_MAX_DIMS)
const gcount       = zeros(UInt,NC_MAX_DIMS)

for (t,ending,aname) in funext
    @eval const $(Symbol(aname)) = zeros($t,1)
end

"""
    nc_char2string(x::Array{UInt8})

Convert a `UInt8` array read from a NetCDF variable of type `NC_CHAR` to a Julia array of Strings
"""
function nc_char2string(x::Array{UInt8})
    y = copy(x)
    y[end,:] .= 0
    dropdims(mapslices(i -> unsafe_string(pointer(i)), y, dims=[1]), dims=1)
end

"""
    nc_string2char(x::Array{String})

Convert a Julia `String` array to an `Array{UInt8}` so that it can be written to a
NetCDF variable of type `NC_CHAR`.
"""
function nc_string2char(s::Array{String})
    maxlen=maximum(length(i) for i in s)
    c = zeros(UInt8,maxlen+1,size(s)...)
    offs=1
    for i=1:length(s)
        copyto!(c,offs,s[i],1)
        offs+=maxlen+1
    end
    c
end

function nc_open(fname::AbstractString,omode::UInt16)
    # Function to open file fname, returns a NetCDF file ID
    nc_open(fname,omode,ida)
    return ida[1]
end

function nc_create(name,mode)
    nc_create(name,mode,ida);
    return ida[1]
end

function nc_inq_dim(id::Integer,idim::Integer)
    # File to inquire dimension idimm returns dimension name and length
    nc_inq_dim(id,idim,namea,lengtha)
    namea[end]=0
    name=unsafe_string(pointer(namea))
    return (name,lengtha[1])
end

function nc_inq_dimid(id::Integer,name::AbstractString)
    # Function to read dimension id for a given function name
    NetCDF.nc_inq_dimid(id,name,dimida)
    return dimida[1]
end

function nc_redef(nc::NcFile)
    if (!nc.in_def_mode)
        nc_redef(nc.ncid)
        nc.in_def_mode=true
    end
end

function nc_enddef(nc::NcFile)
    if nc.in_def_mode
        nc_enddef(nc.ncid)
        nc.in_def_mode=false
    end
end

function nc_inq(id::Integer)
    # Inquire NetCDF file, return number of dims, number of variables, number of global attributes and number of unlimited dimensions
    nc_inq(id,ndima,nvara,ngatta,nunlimdimida)
    return (ndima[1],nvara[1],ngatta[1],nunlimdimida[1])
end

function nc_inq_attname(ncid::Integer, varid::Integer, attnum::Integer)
    # Get attribute name from attribute number
    nc_inq_attname(ncid,varid,attnum,namea)
    namea[end] = 0
    name = unsafe_string(pointer(namea))
    return name
end

function nc_get_att(ncid::Integer, varid::Integer, att::Union{Integer,AbstractString})
    #Reads attribute name, type and number of values
    name = typeof(att)<:Integer ? nc_inq_attname(ncid,varid,att) : string(att)
    nc_inq_att(ncid,varid,name,typea,nvals)
    text = nc_get_att(ncid,varid,string(name),typea[1],nvals[1])
    return name, text
end

#Define methods for writing attributes
nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::Array{Int8})    = nc_put_att_schar(ncid,varid,name,NC_BYTE,length(val),val)
nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::Array{Int16})   = nc_put_att_short(ncid,varid,name,NC_SHORT,length(val),val)
nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::Array{Int32})   = nc_put_att_int(ncid,varid,name,NC_INT,length(val),val)
nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::Array{Int64})   = nc_put_att_long(ncid,varid,name,NC_INT64,length(val),val)
nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::Array{Float32}) = nc_put_att_float(ncid,varid,name,NC_FLOAT,length(val),val)
nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::Array{Float64}) = nc_put_att_double(ncid,varid,name,NC_DOUBLE,length(val),val)
#nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::AbstractString)         = nc_put_att_text(ncid,varid,name,length(val)+1,val)
nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::Array{Any})     = error("Writing attribute array of type Any is not possible")

nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::Int8) = begin int8a[1] = val; nc_put_att(ncid,varid,name,int8a) end
nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::Int16) = begin int16a[1] = val; nc_put_att(ncid,varid,name,int16a) end
nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::Int32) = begin int32a[1] = val; nc_put_att(ncid,varid,name,int32a) end
nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::Int64) = begin int64a[1] = val; nc_put_att(ncid,varid,name,int64a) end
nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::Float32) = begin float32a[1] = val; nc_put_att(ncid,varid,name,float32a) end
nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::Float64) = begin float64a[1] = val; nc_put_att(ncid,varid,name,float64a) end
nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val) = error("Writing attribute of type $(typeof(val)) is currently not possible.")

function nc_put_att(ncid::Integer,varid::Integer,name::AbstractString,val::AbstractString)
    val = string(val)
    len = sizeof(val)
    nc_put_att_text(ncid,varid,name,len,val)
end

function nc_get_att(ncid::Integer,varid::Integer,name::AbstractString,attype::Integer,attlen::Integer)
    if attype == NC_CHAR
        valsa = zeros(UInt8,attlen+1)
        nc_get_att_text(ncid,varid,name,valsa)
        valsa[end] = 0
        return unsafe_string(pointer(valsa))
    else
        valsa = Array{nctype2jltype[attype]}(undef,attlen)
        return nc_get_att!(ncid,varid,name,valsa)
    end
end

nc_get_att!(ncid::Integer,varid::Integer,name::AbstractString,valsa::Array{UInt8}) = begin nc_get_att_ubyte(ncid,varid,name,valsa);valsa end
nc_get_att!(ncid::Integer,varid::Integer,name::AbstractString,valsa::Array{Int8})  = begin nc_get_att_schar(ncid,varid,name,valsa); valsa end
nc_get_att!(ncid::Integer,varid::Integer,name::AbstractString,valsa::Array{Int16})  = begin nc_get_att_short(ncid,varid,name,valsa); valsa end
nc_get_att!(ncid::Integer,varid::Integer,name::AbstractString,valsa::Array{Int32})  = begin nc_get_att_int(ncid,varid,name,valsa); valsa end
nc_get_att!(ncid::Integer,varid::Integer,name::AbstractString,valsa::Array{Int64})  = begin nc_get_att_long(ncid,varid,name,valsa); valsa end
nc_get_att!(ncid::Integer,varid::Integer,name::AbstractString,valsa::Array{Float32})  = begin nc_get_att_float(ncid,varid,name,valsa); valsa end
nc_get_att!(ncid::Integer,varid::Integer,name::AbstractString,valsa::Array{Float64})  = begin nc_get_att_double(ncid,varid,name,valsa); valsa end

function nc_get_att!(ncid::Integer,varid::Integer,name::AbstractString,valsa::Array{T,N}) where {T<:AbstractString,N}
  #Assert that length of the attribute matches length of output
  nc_inq_attlen(ncid, varid, name, lengtha)
  @assert lengtha[1]==length(valsa)
  valsa_c = Array{Ptr{UInt8}}(length(valsa))
  nc_get_att_string(ncid,varid,name,valsa_c)
  for i = 1:length(valsa)
    valsa[i] = unsafe_string(valsa_c[i])
  end
  nc_free_string(length(valsa_c),valsa_c)
  valsa
end

function nc_inq_var(nc::NcFile,varid::Integer)
    # Inquire variables in the file
    nc_inq_var(nc.ncid,varid,namea,typea,ndima,dimida,natta)
    dimids = ndima[1] > 0 ? dimida[1:ndima[1]] : Int32[]
    namea[end] = 0
    name = unsafe_string(pointer(namea))
    #Find out chunks
    nc_inq_var_chunking(nc.ncid,varid,storagep,chunk_sizea)
    chunksize=storagep[1] == NC_CONTIGUOUS ? ntuple(i->0,ndima[1]) : ntuple(i->chunk_sizea[i],ndima[1])
    return name, typea[1], dimids, natta[1], ndima[1], isdimvar(nc,name), chunksize
end

#Test if a variable name is also a dimension name
isdimvar(v::NcVar) = v.name == v.dim[1].name ? true : false

function isdimvar(nc::NcFile,name::AbstractString)
    for n in nc.dim
        if n[2].name==name
            return true
        end
    end
    return false
end

function getdimnamebyid(nc::NcFile,dimid::Integer)
    da = ""
    for d in nc.dim
        da = d[2].dimid == dimid ? d[2].name : da
    end
    return da
end

function getatts_all(ncid::Integer, varid::Integer, natts::Integer)
    atts = Dict{Any,Any}()
    for attnum = 0:natts-1
        attname, attval = nc_get_att(ncid, varid, attnum)
        if ((length(attval)==1) && !(typeof(attval)<:AbstractString))
            attval = attval[1]
        end
        atts[attname] = attval
    end
    NC_VERBOSE && println(atts)
    return atts
end

function getatts(ncid::Integer, varid::Integer, attnames::Union{AbstractString,AbstractArray})
    atts = Dict{Any,Any}()
    for aattname = (typeof(attnames)<:AbstractArray ? attnames : [attnames])
        attname, attval = nc_get_att(ncid, varid, aattname)
        if ((length(attval)==1) && !(typeof(attval)<:AbstractString))
            attval = attval[1]
        end
        atts[attname] = attval
    end
    NC_VERBOSE && println(atts)
    return atts
end

function readdimvars(nc::NcFile)
    for v in nc.vars
        if isdimvar(v[2])
            v[2].dim[1].vals = readvar(nc,v[1])
            d[2].atts = v[2].atts
        end
    end
end

function checkboundsNC(v)
    nd = ndims(v)
    for i = 1:nd
        ci = nd + 1 - i
        gstart[ci] < 0 && error("Start index must not be smaller than 1")
        if gstart[ci] + gcount[ci] > v.dim[i].dimlen
            if v.dim[i].unlim
                #Reset length of unlimited dimension
                v.dim[i].dimlen=gstart[ci] + gcount[ci]
            else
                error("Start + Count exceeds dimension length in dimension $(v.dim[i].name)")
            end
        end
        nothing
    end
end

function preparestartcount(start, count, v::NcVar)

    length(start) == v.ndim || error("Length of start ($(length(start))) must equal the number of variable dimensions ($(v.ndim))")
    length(count) == v.ndim || error("Length of start ($(length(count))) must equal the number of variable dimensions ($(v.ndim))")

    p = one(eltype(gcount))
    nd = length(start)

    for i=1:nd
        ci = nd + 1 - i
        gstart[ci] = start[i] - 1
        gcount[ci] = count[i] < 0 ? v.dim[i].dimlen - gstart[ci] : count[i]

        p *= gcount[ci]
    end

    checkboundsNC(v)

    return p
end

function _readdimvars(nc::NcFile)
    for d in nc.dim
        for v in nc.vars
            if d[2].name == v[2].name
                NC_VERBOSE && println(d[2].name," ",v[2].name)
                d[2].vals = readvar(nc,v[2].name)
                d[2].atts = v[2].atts
            end
        end
    end
end

ischunked(v::NcVar) = v.chunksize[1] > 0

defaultstart(v::NcVar) = ones(Int, v.ndim)
defaultcount(v::NcVar) = Int[i for i in size(v)]

function parsedimargs(dim)
  isempty(dim) && return NcDim[]
    idim = 0
    dimlen = nothing
    dimvals = nothing
    dimatts = nothing
    name = nothing
    d = NcDim[]
    for a in dim
        NC_VERBOSE && println(a,idim)
        if isa(a, AbstractString)
            #Assume a name is given
            #first create an NcDim object from the last dim
            if (name != nothing)
                push!(d, finalizedim(dimlen,dimvals,dimatts,name))
                idim = idim+1
                dimlen = nothing
                dimvals = nothing
                dimatts = nothing
                name = nothing
            end
            name = a
        elseif isa(a, Integer)
            #Assume a dimension length is given
            if a > 0
                dimlen = a
            else
                #Assume unlimited dimension
                dimlen = 0
            end
        elseif isa(a, AbstractFloat) && isinf(a)
            #Generate unlimited dimension
            dimlen = 0
        elseif isa(a, AbstractArray)
            #Assume dimension values are given
            if dimvals == nothing
                dimvals = a
                dimlen = length(dimvals)
            else
                error("Dimension values of $name defined more than once")
            end
        elseif isa(a, Dict)
            #Assume attributes are given
            dimatts = dimatts == nothing ? a : error("Dimension attributes of $name defined more than once")
        else
            error("Could not parse argument $a in nccreate")
        end
    end
    push!(d, finalizedim(dimlen,dimvals,dimatts,name))
    return d
end

function finalizedim(dimlen, dimvals, dimatts, name)
    if (dimlen == nothing) && (dimvals == nothing)
        dimlen = 1
    end
    if (dimlen != nothing) && (dimvals == nothing)
        dimvals = Array{Float64}(undef,0)
    end
    if dimatts == nothing
        dimatts = Dict("missval"=>-9999)
    end
    isunlimited = dimlen == 0 ? true : false
    return NcDim(name, dimlen, atts=dimatts, values=dimvals, unlimited=isunlimited)
end


text:
using CUDA
using CUDA.CUSPARSE
using ExaPF
using KernelAbstractions
using LinearAlgebra
using Random
using SparseArrays
using Test
using TimerOutputs
using Krylov

const LS = ExaPF.LinearSolvers

@testset "Iterative linear solvers with custom block Jacobi" begin
    # Square and preconditioned problems.
    function square_preconditioned(n :: Int=10)
        A   = ones(n, n) + (n-1) * Matrix(I, n, n)
        b   = 10.0 * [1:n;]
        M⁻¹ = 1/n * Matrix(I, n, n)
        return A, b, M⁻¹
    end
    A, b, N = square_preconditioned(100)
    x = similar(b); r = similar(b)
    spA = sparse(A)
    nblocks = 2
    precond = LS.BlockJacobiPreconditioner(spA, nblocks, CPU())
    LS.update(spA, precond)
    @testset "BICGSTAB" begin
        LS.ldiv!(LS.BICGSTAB(precond), x, spA, b)
        r = b - spA * x
        resid = norm(r) / norm(b)
        @test(resid ≤ 1e-6)
    end
    # Embed preconditioner in linear solvers
    @testset "($LinSolver)" for LinSolver in ExaPF.list_solvers(CPU())
        algo = LinSolver(precond)
        LS.ldiv!(algo, x, spA, b)
        r = b - spA * x
        resid = norm(r) / norm(b)
        @test(resid ≤ 1e-6)
    end
end

@testset "Wrapping of iterative solvers" begin
    nblocks = 2
    n, m = 32, 32

    # Add a diagonal term for conditionning
    A = randn(n, m) + 15I
    x♯ = randn(m)
    b = A * x♯

    # Be careful: all algorithms work with sparse matrix
    A = sparse(A)
    # Init iterators for KernelAbstractions
    if has_cuda_gpu()
        iterators = zip([CPU(), CUDADevice()],
                        [Array{Float64}, CuVector{Float64}],
                        [SparseMatrixCSC, CuSparseMatrixCSR])
    else
        iterators = zip([CPU()], [Array{Float64}], [SparseMatrixCSC])
    end
    @testset "Run iterative solvers on device $(device)" for (device, V, SM) in iterators
        As = SM(A)
        bs = convert(V, b)
        x0 = convert(V, zeros(m))
        xs♯ = convert(V, x♯)

        # TODO: currently Preconditioner takes as input only sparse matrix
        # defined on the main memory.
        precond = LS.BlockJacobiPreconditioner(A, nblocks, device)

        # First test the custom implementation of BICGSTAB
        @testset "BICGSTAB" begin
            # Need to update preconditioner before resolution
            LS.update(As, precond)
            fill!(x0, 0.0)
            n_iters = LS.ldiv!(LS.BICGSTAB(precond), x0, As, bs)
            @test n_iters <= m
            @test x0 ≈ xs♯ atol=1e-6
        end
        @testset "Interface for iterative algorithm ($LinSolver)" for LinSolver in ExaPF.list_solvers(device)
            algo = LinSolver(precond)
            fill!(x0, 0.0)
            n_iters = LS.ldiv!(algo, x0, As, bs)
            @test n_iters <= m
            @test x0 ≈ xs♯ atol=1e-6
        end
    end
end


text:
# Convert the DFT data from https://github.com/FitSNAP/FitSNAP/tree/master/examples/Ta_Linear_JCP2014/JSON
# to the current format required by PotentialLearning.jl.
# The current format will probably change in the coming months.

using JSON

path = "./JSON/"
for d in readdir(path)
    sub_path = joinpath(path, d, "DATA")
    rm(sub_path, recursive=true)
    dirs = readdir(joinpath(path, d))
    mkdir(sub_path)
    for (i, f) in enumerate(dirs)
        global dft_data = Dict()
        open(joinpath(path, d, f) , "r") do file
            @show joinpath(path, d, f)
            readline(file)
            dft_data = JSON.parse(readline(file))
        end
        ssub_path = joinpath(sub_path, "$i")
        mkdir(ssub_path)

        open(joinpath(ssub_path, "ENERGY"), "w") do file
            write(file, "$(dft_data["Dataset"]["Data"][1]["Energy"])\n")
        end

        open(joinpath(ssub_path, "FORCES"), "w") do file
            for force in dft_data["Dataset"]["Data"][1]["Forces"]
                write(file, "$(force[1]) $(force[2]) $(force[3])\n")
            end
        end

        open(joinpath(ssub_path, "DATA") , "w") do file
            write(file, "LAMMPS DATA file\n")
            write(file, "\n")
            write(file, "$(dft_data["Dataset"]["Data"][1]["NumAtoms"]) atoms\n")
            write(file, "0 bonds\n")
            write(file, "0 angles\n")
            write(file, "0 dihedrals\n")
            write(file, "0 impropers\n")
            write(file, "\n")
            write(file, "1 atom types\n")
            write(file, "0 bond types\n")
            write(file, "0 angle types\n")
            write(file, "\n")
            write(file, "0.0 12.990503 xlo xhi\n")
            write(file, "0.0 11.250148 ylo yhi\n")
            write(file, "0.0 15.845291 zlo zhi\n")
            write(file, "\n")
            write(file, "Masses\n")
            write(file, "\n")
            write(file, "1 73\n")
            write(file, "\n")
            write(file, "Atoms\n")
            write(file, "\n")
            positions = dft_data["Dataset"]["Data"][1]["Positions"]
            for (j, pos) in enumerate(positions)
                write(file, "$j 1 $(pos[1]) $(pos[2]) $(pos[3]) 0 0 0\n")
            end
        end
    end
end


text:
########################################################################################################
## Details
# This file implements the MAP estimation proceedure for UQ inference problems.
#       Note, we convert the named tuple p.x into a vector for use in GalacticOptim.
#       For the conversion, we use TransformVariables.
#########################################################################################################

function MAP(pdist :: ArbitraryDistribution; solver = LBFGS(), verbose = false)

    # Define a Hamiltonian system
    x0 = inverse(pdist.t, pdist.x)
    n = length(x0)
    @model function demo()
        x ~ pdist.prior
        Turing.@addlogprob! pdist.distribution(x)
    end

    model = demo()
    map_estimate = optimize(model, Turing.MAP(), x0, solver)
    pdist.x = transform(pdist.t, map_estimate.values.array)
    if verbose
        println("Map Estimate ")
        show(stdout, "text/plain", map_estimate)
    end

    return pdist

end


text:
function add_kepler_contribution(obj::FieldTimeDerivativesEquations, arg0::RealFieldElement)
    return jcall(obj, "addKeplerContribution", void, (RealFieldElement,), arg0)
end

function add_mass_derivative(obj::FieldTimeDerivativesEquations, arg0::RealFieldElement)
    return jcall(obj, "addMassDerivative", void, (RealFieldElement,), arg0)
end

function add_non_keplerian_acceleration(obj::FieldTimeDerivativesEquations, arg0::FieldVector3D)
    return jcall(obj, "addNonKeplerianAcceleration", void, (FieldVector3D,), arg0)
end


text:
using GymSpace: Box, Discrete

"""
Copyright 2013, RLPy http://acl.mit.edu/RLPy
Credit: Alborz Geramifard, Robert H. Klein, Christoph Dann,
               William Dabney, Jonathan P. How
License: "BSD 3-Clause"
Author: "Christoph Dann <cdann@cdann.de>"
"""
mutable struct AcrobotEnv <: AbstractEnv
    viewer
    observation_space::Box
    action_space::Discrete
    state
    dt::Float32
    LINK_LENGTH_1::Float32  # [m]
    LINK_LENGTH_2::Float32  # [m]
    LINK_MASS_1::Float32  # [kg] mass of link 1
    LINK_MASS_2::Float32  # [kg] mass of link 2
    LINK_COM_POS_1::Float32  # [m] position of the center of mass of link 1
    LINK_COM_POS_2::Float32  # [m] position of the center of mass of link 2
    LINK_MOI::Float32  # moments of intertia for both links
    MAX_VEL_1::Float32
    MAX_VEL_2::Float32
    AVAIL_TORQUE::Array{Float32, N} where N
    torque_noise_max::Float32
    book_or_nips::String
end

include("vis/acrobot.jl")

function AcrobotEnv()
    dt = 2f-1
    LINK_LENGTH_1 = 1f0
    LINK_LENGTH_2 = 1f0
    LINK_MASS_1 = 1f0
    LINK_MASS_2 = 1f0
    LINK_COM_POS_1 = 5f-1
    LINK_COM_POS_2 = 5f-1
    LINK_MOI = 1f0

    MAX_VEL_1 = Float32(4 * pi)
    MAX_VEL_2 = Float32(9 * pi)

    AVAIL_TORQUE = [-1f0, 0f0, 1f0]
    torque_noise_max = 0f0

    # use dynamics equations from the nips paper or the book
    book_or_nips = "book"
    # action_arrow = nothing
    # domain_fig = nothing
    # actions_num = 3

    viewer = nothing
    high = [1f0, 1f0, 1f0, 1f0, MAX_VEL_1, MAX_VEL_2] #//
    low = -high
    observation_space = Box(low, high, Float32)
    action_space = Discrete(3)
    state = nothing

    AcrobotEnv(viewer, observation_space, action_space, state, dt, LINK_LENGTH_1,
               LINK_LENGTH_2, LINK_MASS_1, LINK_MASS_2, LINK_COM_POS_1, LINK_COM_POS_2,
               LINK_MOI, MAX_VEL_1, MAX_VEL_2, AVAIL_TORQUE, torque_noise_max, book_or_nips)
end

function reset!(env::AcrobotEnv)
    env.state = param(rand(Float32, 4) * 2f-1 .- 1f-1)
    return _get_ob(env)
end

#=
function custom_reset!(env::AcrobotEnv, args=[-0.04616826f0, -0.04653378f0,  0.07321809f0, -0.04461966f0])
    env.state = param(args)
    return _get_ob(env)
end
=#
get_torque(env::AcrobotEnv, action) = env.AVAIL_TORQUE[action]

function step!(env::AcrobotEnv, action)
    @assert action ∈ env.action_space "$(action) is not in the action space. Enter a valid action."
    s = env.state
    torque = get_torque(env, action)

    # Add noise to the force action
    env.torque_noise_max > 0 && (torque += env.torque_noise_max * (2rand(Float32) - 1))

    # Now, augment the state with our force action so it can be passed to _dsdt
    s_augmented = vcat(s, torque)

    ns = rk4(env, _dsdt, s_augmented, [0, env.dt])
    # only care about final timestep of integration returned by integrator
    ns = ns[end, :]

    ns = vcat(wrap(ns[1:2], -π, π), bound(ns[3:4], -env.MAX_VEL_1, env.MAX_VEL_2)) # omit action
    env.state = ns
    terminal = _terminal(env)
    reward = !all(terminal) ? -1f0 : 0f0

    return _get_ob(env), reward, terminal, Dict()
end

function _get_ob(env::AcrobotEnv)
    s = env.state
    return vcat(cos.(s[1:1]), sin.(s[1:1]), cos.(s[2:2]), sin.(s[2:2]), s[3:4])
end

function _terminal(env::AcrobotEnv)
    s = env.state
    return cos(s[1]) - cos(s[2] + s[1]) > 1f0
end

function _dsdt(env::AcrobotEnv, s_augmented, t)
    m1 = env.LINK_MASS_1
    m2 = env.LINK_MASS_2
    l1 = env.LINK_LENGTH_1
    lc1 = env.LINK_COM_POS_1
    lc2 = env.LINK_COM_POS_2
    I1 = env.LINK_MOI
    I2 = env.LINK_MOI
    g = 98f-1
    a = s_augmented[end:end]
    s = s_augmented[1:end-1]
    theta1 = s[1:1]
    theta2 = s[2:2]
    dtheta1 = s[3:3]
    dtheta2 = s[4:4]

    d1 = m1 * lc1^2 .+ m2 * (l1^2 .+ lc2^2 .+ 2l1 * lc2 * cos.(theta2)) .+ I1 .+ I2

    d2 = m2 * (lc2^2 .+ l1 * lc2 * cos.(theta2)) .+ I2

    phi2 = m2 * lc2 * g * cos.(theta1 .+ theta2 .- pi/2f0)

    phi1 = - m2 * l1 * lc2 * dtheta2.^2 .* sin.(theta2) .-
            2 * m2 * l1 * lc2 * dtheta2 .* dtheta1 .* sin.(theta2) .+
            (m1 * lc1 + m2 * l1) * g * cos.(theta1 .- Float32(pi) / 2) .+ phi2

    if env.book_or_nips == "nips"
        # the following line is consistent with the description in the paper
        ddtheta2 = (a .+ d2./d1 .* phi1 .- phi2) ./
                    (m2 .* lc2^2 .+ I2 .- d2.^2 ./d1)
    else
        # the following line is consistent with the java impelemtation and the
        # book
        ddtheta2 = (a .+ d2 ./ d1 .* phi1 .- m2 .* l1 .* lc2 .* dtheta1.^2 .* sin.(theta2) .-phi2) ./
                    (m2 .* lc2^2 .+ I2 .- d2.^2 ./ d1)
    end
    ddtheta1 = -(d2 .* ddtheta2 .+ phi1) ./ d1

    return vcat(dtheta1, dtheta2, ddtheta1, ddtheta2, 0f1)
end

#=
x : value to be wrapped\n
m : minimum possible value in range\n
M : maximum possible value in range\n
Wraps ``x`` so m <= x <= M; but unlike ``bound()`` which truncates,
``wrap()`` wraps x around the coordinate system defined by m, M.\n
For example, m = -180, M = 180 (degrees), x = 360 --> returns 0.
=#
function wrap(x, m, M)
    diff = M - m
    while all(x .> M)
        x .-= diff
    end
    while all(x .< m)
        x .+= diff
    end
    return x
end

# returns m <= x <= M
bound(x, m, M) = min.(max.(x, m), M)

#=
y0 : initial state vector
t : sample times
derivs : returns the derivative of the system and has the
         signature ``dy = derivs(yi, ti)``
args : additional arguments passed to the derivative function
kwargs : additional keyword arguments passed to the derivative function
=#
function rk4(env::AcrobotEnv, derivs, y0, t, args...; kwargs...)
    yout = param(zeros(Float32, length(t), length(y0)))
    yout.data[1, :] = y0.data

    for i=1:(length(t)-1)
        thist = t[i]
        dt = t[i+1] - thist
        dt2 = dt / 2f0
        y0 = yout[i, :]

        k1 = derivs(env, y0, thist, args...; kwargs...)
        k2 = derivs(env, y0 .+ dt2 * k1, thist .+ dt2, args...; kwargs...)
        k3 = derivs(env, y0 .+ dt2 * k2, thist .+ dt2, args...; kwargs...)
        k4 = derivs(env, y0 .+ dt * k3, thist .+ dt, args...; kwargs...)

        yout.data[i+1, :] = (y0 .+ dt / 6f0 * (k1 .+ 2 * k2 .+ 2 * k3 .+ k4)).data

    end
    return yout
end

show(io::IO, env::AcrobotEnv) = print(io, "AcrobotEnv")


text:
import Base: axpy!, scale!, ctranspose, dot, size
import ArrayViews: view

export AbstractSketch, IdentitySketch, SymmetricSketch, AsymmetricSketch,
	additive_update!, reconstruct, dot

abstract AbstractSketch

# Identity

type IdentitySketch<:AbstractSketch
	X
end

function additive_update!(s::IdentitySketch, tilde_X, eta=1)
	axpy!(eta, tilde_X, s.X)
end

function cgd_update!(s::IdentitySketch, tilde_X, eta)
	scale!(s.X, 1-eta)
	axpy!(eta, tilde_X, s.X)
end

function reconstruct(s::IdentitySketch)
	return s.X
end

# Symmetric

type SymmetricSketch<:AbstractSketch
	r::Int
	Omega::Array{Float64,2}
	Psi::Array{Float64,2}
	Y::Array{Float64,2}
	W::Array{Float64,2}
end

function SymmetricSketch(n::Int, r::Int)
	k = 2r + 1
	l = 4r + 3 # = 2k+1
	Omega = randn(n,k)
	Psi   = randn(n,l)
	Y     = zeros(n,k)
	W     = zeros(n,l)
	return SymmetricSketch(r, Omega, Psi, Y, W)
end

function additive_update!(s::SymmetricSketch, vvt, eta=1)
	# Y = Y + eta*v*v'*Omega
	axpy!(eta, vvt*s.Omega, s.Y)
	# W = W + eta*v*v'*Psi
	axpy!(eta, vvt*s.Psi, s.W)
end

function cgd_update!(s::SymmetricSketch, vvt, eta)
	# Y = (1-eta)Y + eta*v*v'*Omega
	scale!(s.Y, 1-eta)
	axpy!(eta, vvt*s.Omega, s.Y)
	# W = (1-eta)W + eta*v*v'*Psi
	scale!(s.W, 1-eta)
	axpy!(eta, vvt*s.Psi, s.W)
end

function reconstruct(s::SymmetricSketch)
	# Q = orth(s.Y)
	Q,_ = qr(s.Y)
	B = (Q's.W) / (Q's.Psi)
	d,V = eigs((B+B')/2, nev=s.r, which=:LM)
	U = Q*V
	d = pos(d)
	return LowRankOperator(U, spdiagm(d)*U') # reconstruction as square matrix is U*diagm(d)*U'
end

# Asymmetric

type AsymmetricSketch<:AbstractSketch
	r::Int
	Omega::Array{Float64,2}
	Psi::Array{Float64,2}
	Y::Array{Float64,2}
	W::Array{Float64,2}
end

function AsymmetricSketch(m::Int, n::Int, r::Int)
	k = 2r + 1
	l = 4r + 3 # == 2k+1
	Omega = randn(n,k)
	Psi   = randn(m,l)
	Y     = zeros(m,k)
	W     = zeros(n,l)
	return AsymmetricSketch(r, Omega, Psi, Y, W)
end

function additive_update!(s::AsymmetricSketch, uvt, eta=1)
	# Y = Y + eta*v*v'*Omega
	axpy!(eta, uvt*s.Omega, s.Y)
	# W = W + eta*v*v'*Psi
	axpy!(eta, uvt'*s.Psi, s.W)
end

function cgd_update!(s::AsymmetricSketch, uvt, eta)
	# Y = (1-eta)Y + eta*v*v'*Omega
	scale!(s.Y, 1-eta)
	axpy!(eta, uvt*s.Omega, s.Y)
	# W = (1-eta)W + eta*v*v'*Psi
	scale!(s.W, 1-eta)
	axpy!(eta, uvt'*s.Psi, s.W)
end

function reconstruct(s::AsymmetricSketch, r::Int=s.r)
	# Q = orth(s.Y)
	k = 2r + 1
	l = 4r + 3
	Q,_ = qr(view(s.Y,:,1:k))
	B = view(s.W,:,1:l) / (Q'*view(s.Psi,:,1:l)) # Q's.Psi is k x l, its pinv is l x k, so B is n x k
	U,s,V = mysvds(B, nsv=r) # U is n x r
	return LowRankOperator(Q*V, spdiagm(s)*U') # reconstruction as square matrix is Q*V*diagm(s)*U'
end
# function reconstruct(s::AsymmetricSketch)
# 	# Q = orth(s.Y)
# 	Q,_ = qr(s.Y)
# 	B = s.W * pinv(Q'*s.Psi) # B is k x l
# 	U,s,V,_ = svds(B, nsv=s.r)
# 	return LowRankOperator(Q*U, spdiagm(s), V) # reconstruction as square matrix is Q*U*diagm(s)*V'
# end

### utilities
pos(a::Number) = max(a,0)
function pos(a::Array)
	for i in length(a)
		a[i] = pos(a[i])
	end
	a
end

# not yet efficient
function dot(M::SparseMatrixCSC{Float64,Int64}, A::Array{Float64,2})
	dot(M.nzval, A[find(M)])
end

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
module DataBench

using PooledArrays, RCall
import Base.ht_keyindex
# using FastGroupBy

export R_bench, createIndexedTable, run_juliadb_bench, run_juliadb_bench_pmap
export gen_string_vec_fixed_len, gen_string_vec_var_len, gen_string_vec_id_fixed_len
export createSynDataFrame

# include("julia/blogpost/indexedtable.jl_bench.jl")
include("julia/r_data_table_group_by/rbench.jl")
include("julia/r_data_table_group_by/DataFrames.jl")

end # module


text:
@testset "others" begin
    @testset "TuringWishart" begin
        dim = 3
        A = Matrix{Float64}(I, dim, dim)
        dW1 = Wishart(dim + 4, A)
        dW2 = TuringWishart(dim + 4, A)
        dW3 = TuringWishart(dW1)
        mean = Distributions.mean
        @testset "$F" for F in (size, rank, mean, meanlogdet, entropy, cov, var)
            @test F(dW1) == F(dW2) == F(dW3)
        end
        @test Matrix(mode(dW1)) == mode(dW2) == mode(dW3)
        xw = rand(dW2)
        @test insupport(dW1, xw)
        @test insupport(dW2, xw)
        @test insupport(dW3, xw)
        @test logpdf(dW1, xw) == logpdf(dW2, xw) == logpdf(dW3, xw)
    end

    @testset "TuringInverseWishart" begin
        dim = 3
        A = Matrix{Float64}(I, dim, dim)
        dIW1 = InverseWishart(dim + 4, A)
        dIW2 = TuringInverseWishart(dim + 4, A)
        dIW3 = TuringInverseWishart(dIW1)
        mean = Distributions.mean
        @testset "$F" for F in (size, rank, mean, cov, var)
            @test F(dIW1) == F(dIW2) == F(dIW3)
        end
        @test Matrix(mode(dIW1)) == mode(dIW2) == mode(dIW3)
        xiw = rand(dIW2)
        @test insupport(dIW1, xiw)
        @test insupport(dIW2, xiw)
        @test insupport(dIW3, xiw)
        @test logpdf(dIW1, xiw) == logpdf(dIW2, xiw) == logpdf(dIW3, xiw)
    end

    @testset "TuringMvNormal" begin
        @testset "$TD" for TD in [TuringDenseMvNormal, TuringDiagMvNormal, TuringScalMvNormal]
            m = rand(3)
            if TD <: TuringDenseMvNormal
                A = rand(3, 3)
                C = A' * A + I
                d1 = TuringMvNormal(m, C)
                d2 = MvNormal(m, C)
            elseif TD <: TuringDiagMvNormal
                C = rand(3)
                d1 = TuringMvNormal(m, C)
                d2 = MvNormal(m, Diagonal(C .^ 2))
            else
                C = rand()
                d1 = TuringMvNormal(m, C)
                d2 = MvNormal(m, C^2 * I)
            end

            @testset "$F" for F in (length, size, mean)
                @test F(d1) == F(d2)
            end
            @test cov(d1) ≈ cov(d2)
            @test var(d1) ≈ var(d2)

            x1 = rand(d1)
            x2 = rand(d1, 3)
            @test isapprox(logpdf(d1, x1), logpdf(d2, x1), rtol = 1e-6)
            @test isapprox(logpdf(d1, x2), logpdf(d2, x2), rtol = 1e-6)
        end
    end

    @testset "TuringMvLogNormal" begin
        @testset "$TD" for TD in [TuringDenseMvNormal, TuringDiagMvNormal, TuringScalMvNormal]
            m = rand(3)
            if TD <: TuringDenseMvNormal
                C = Matrix{Float64}(I, 3, 3)
                d1 = TuringMvLogNormal(TuringMvNormal(m, C))
                d2 = MvLogNormal(MvNormal(m, C))
            elseif TD <: TuringDiagMvNormal
                C = ones(3)
                d1 = TuringMvLogNormal(TuringMvNormal(m, C))
                d2 = MvLogNormal(MvNormal(m, Diagonal(C .^ 2)))
            else
                C = 1.0
                d1 = TuringMvLogNormal(TuringMvNormal(m, C))
                d2 = MvLogNormal(MvNormal(m, C^2 * I))
            end

            @test length(d1) == length(d2)

            x1 = rand(d1)
            x2 = rand(d1, 3)
            @test logpdf(d1, x1) ≈ logpdf(d2, x1) rtol=1e-6
            @test logpdf(d1, x2) ≈ logpdf(d2, x2) rtol=1e-6

            x2[:, 1] .= -1
            @test isinf(logpdf(d1, x2)[1])
            @test isinf(logpdf(d2, x2)[1])
        end
    end

    @testset "TuringUniform" begin
        @test logpdf(TuringUniform(), 0.5) == 0
    end

    @testset "TuringPoissonBinomial" begin
        d1 = TuringPoissonBinomial([0.5, 0.5])
        d2 = PoissonBinomial([0.5, 0.5])
        @test quantile(d1, 0.5) == quantile(d2, 0.5)
        @test minimum(d1) == minimum(d2)
    end

    @testset "Inverse of pi" begin
        @test 1/pi == inv(pi)
    end

    @testset "Cholesky" begin
        A = rand(3, 3)'; A = A + A' + 3I;
        C = cholesky(A; check = true)
        factors, info = DistributionsAD.turing_chol(A, true)
        @test factors == C.factors
        @test info == C.info
    end

    @testset "zygote_ldiv" begin
        A = to_posdef(rand(3, 3))
        B = copy(A)
        @test DistributionsAD.zygote_ldiv(A, B) == A \ B
    end

    @testset "Entropy" begin
        sigma = exp(randn())
        d1 = TuringScalMvNormal(randn(10), sigma)
        d2 = MvNormal(randn(10), sigma^2 * I)
        @test entropy(d1) ≈ entropy(d2) rtol=1e-6

        sigmas = exp.(randn(10))
        d1 = TuringDiagMvNormal(randn(10), sigmas)
        d2 = MvNormal(randn(10), Diagonal(sigmas .^ 2))
        @test entropy(d1) ≈ entropy(d2) rtol=1e-6

        A = randn(10)
        C = A * A' + I
        d1 = TuringDenseMvNormal(randn(10), C)
        d2 = MvNormal(randn(10), C)
        @test entropy(d1) ≈ entropy(d2) rtol=1e-6
    end

    @testset "Params" begin
        m = rand(10)
        sigmas = randexp(10)

        d = TuringDiagMvNormal(m, sigmas)
        @test params(d) == (m, sigmas)

        d = TuringScalMvNormal(m, sigmas[1])
        @test params(d) == (m, sigmas[1])
    end

    @testset "adapt_randn" begin
        rng = MersenneTwister()
        for T in (Float32, Float64)
            test_adapt_randn(rng, rand(rng, T, 50), T, 10, 30)
        end
    end

    @testset "TuringDirichlet" begin
        dim = 3
        n = 4
        for alpha in (2, rand())
            d1 = TuringDirichlet(dim, alpha)
            d2 = Dirichlet(dim, alpha)
            d3 = TuringDirichlet(d2)
            @test d1.alpha == d2.alpha == d3.alpha
            @test d1.alpha0 == d2.alpha0 == d3.alpha0
            @test d1.lmnB == d2.lmnB == d3.lmnB

            s1 = rand(d1)
            @test s1 isa Vector{Float64}
            @test length(s1) == dim

            s2 = rand(d1, n)
            @test s2 isa Matrix{Float64}
            @test size(s2) == (dim, n)
        end

        for alpha in (ones(Int, dim), rand(dim))
            d1 = TuringDirichlet(alpha)
            d2 = Dirichlet(alpha)
            d3 = TuringDirichlet(d2)
            @test d1.alpha == d2.alpha == d3.alpha
            @test d1.alpha0 == d2.alpha0 == d3.alpha0
            @test d1.lmnB == d2.lmnB == d3.lmnB

            s1 = rand(d1)
            @test s1 isa Vector{Float64}
            @test length(s1) == dim

            s2 = rand(d1, n)
            @test s2 isa Matrix{Float64}
            @test size(s2) == (dim, n)
        end

        # https://github.com/TuringLang/DistributionsAD.jl/issues/158
        let
            d = TuringDirichlet(rand(2))
            z = rand(d)
            logpdf_z = logpdf(d, z)
            pdf_z = pdf(d, z)

            for x in ([0.5, 0.8], [-0.5, 1.5])
                @test logpdf(d, x) == -Inf
                @test iszero(pdf(d, x))

                xmat = hcat(x, x)
                @test all(==(-Inf), logpdf(d, xmat))
                @test all(iszero, pdf(d, xmat))

                xzmat = hcat(x, z)
                @test logpdf(d, xzmat) == [-Inf, logpdf_z]
                @test pdf(d, xzmat) == [0, pdf_z]
            end
        end
    end
end


text:
for SpatialVector in (:Twist, :SpatialAcceleration, :Momentum, :Wrench)
    @eval begin
        # Basics
        Base.eltype(::Type{$SpatialVector{T}}) where {T} = T
        angular(v::$SpatialVector) = v.angular
        linear(v::$SpatialVector) = v.linear

        # Construct/convert given another spatial vector
        $SpatialVector{T}(v::$SpatialVector{T}) where {T} = v
        Base.convert(::Type{V}, v::$SpatialVector) where {V<:$SpatialVector} = V(v)

        # Construct/convert to SVector
        StaticArrays.SVector{6, T}(v::$SpatialVector) where {T} = convert(SVector{6, T}, [angular(v); linear(v)])
        StaticArrays.SVector{6}(v::$SpatialVector{X}) where {X} = SVector{6, X}(v)
        StaticArrays.SVector(v::$SpatialVector) = SVector{6}(v)
        StaticArrays.SArray(v::$SpatialVector) = SVector(v)
        Base.convert(::Type{SA}, v::$SpatialVector) where {SA <: SArray} = SA(v)

        function Base.Array(v::$SpatialVector)
            Base.depwarn("Array(v) has been deprecated. Please use SVector(v) instead.", :Array)
            SVector(v)
        end
    end
end

for SpatialMatrix in (:GeometricJacobian, :MomentumMatrix, :WrenchMatrix)
    @eval begin
        # Basics
        Base.eltype(::Type{$SpatialMatrix{A}}) where {T, A<:AbstractMatrix{T}} = T
        Base.size(m::$SpatialMatrix) = (6, size(angular(m), 2))
        Base.size(m::$SpatialMatrix, d) = size(m)[d]
        Base.transpose(m::$SpatialMatrix) = Transpose(m)
        angular(m::$SpatialMatrix) = m.angular
        linear(m::$SpatialMatrix) = m.linear

        # Construct/convert given another spatial matrix
        $SpatialMatrix(m::$SpatialMatrix{A}) where {A} = SpatialMatrix{A}(m)
        Base.convert(::Type{$SpatialMatrix{A}}, m::$SpatialMatrix) where {A} = $SpatialMatrix{A}(m)
        Base.convert(::Type{$SpatialMatrix{A}}, m::$SpatialMatrix{A}) where {A} = m

        # Construct/convert to Matrix
        (::Type{A})(m::$SpatialMatrix) where {A<:Array} = convert(A, [angular(m); linear(m)])
        Base.convert(::Type{A}, m::$SpatialMatrix) where {A<:Array} = A(m)
    end
end


text:
using Oceananigans.Advection
using Oceananigans.Coriolis
using Oceananigans.Operators

"""
Compute the tendency for the x-directional transport, uh
"""
@inline function uh_solution_tendency(i, j, k, grid,
                                      advection,
                                      coriolis,
                                      solution,
                                      tracers,
                                      diffusivities,
                                      forcings,
                                      clock)

    return ( - 0.0 )
#    return ( - ∂xᶠᵃᵃ(i, j, k, grid, solution.h) )
end

"""
Compute the tendency for the y-directional transport, vh.
"""
@inline function vh_solution_tendency(i, j, k, grid,
                                      advection,
                                      coriolis,
                                      solution,
                                      tracers,
                                      diffusivities,
                                      forcings,
                                      clock)

    return ( -0.0 )
#    return ( - ∂yᵃᶠᵃ(i, j, k, grid, solution.h) )

end

"""
Compute the tendency for the height, h.
"""
@inline function h_solution_tendency(i, j, k, grid,
                                     advection,
                                     coriolis,
                                     solution,
                                     tracers,
                                     diffusivities,
                                     forcings,
                                     clock) where tracer_index

    return ( - 0.0  )
#    return ( - ∂xᶠᵃᵃ(i, j, k, grid, solution.h) )

end

@inline function tracer_tendency(i, j, k, grid,
                                 val_tracer_index::Val{tracer_index},
                                 advection,
                                 solution,
                                 tracers,
                                 diffusivities,
                                 forcings,
                                 clock) where tracer_index

    @inbounds c = tracers[tracer_index]

    # velocities should be replaced with transport, as a vector
    return ( 0.0 )
#    return ( - div_Uc(i, j, k, grid, advection, solution, c) )

end


text:
#cfg title = "Creating a MuJoCo Environment"
#cfg weight = 11
#cfg active = true
#cfg deps = ["humanoid.xml"]

# ## Overview

# Using LyceumMuJoCo, we will create the environment for a Humanoid "get-up" task
# that mostly relies on the defaults of LyceumBase and LyceumMuJoCo to
# propagate state, action, and observation data. We will have to implement reward
# and evaluation functions, of course, along with a few other parts of the
# `AbstractEnvironment` interface.

# We then solve the "get-up" task using a Model-Predictive Control method called
# ["Model Predictive Path Integral Control"](https://ieeexplore.ieee.org/iel7/7478842/7487087/07487277.pdf)
# or MPPI, walking through how to log experiment data and plot the results.

# ## The Code

# First we grab our dependencies of the Lyceum ecosystem and other helpful packages.
using LinearAlgebra, Random, Statistics
using Plots, UnicodePlots, JLSO
using LyceumBase, LyceumBase.Tools, LyceumAI, LyceumMuJoCo, MuJoCo, UniversalLogger, Shapes

#md # ### Humanoid Type

# This struct is our primary entry into the environment API. Environments utilizing
# the MuJoCo simulator through LyceumMuJoCo should subtype
# `AbstractMuJoCoEnvironment <: AbstractEnvironment`. As you can see, this simple
# example only wraps around the underlying simulator (the `sim::MJSim` field of `Humanoid`,
# referred to hereafter as just `sim`). The functions of the LyceumBase API will then
# dispatch on this struct through Julia's [multiple dispatch](https://en.wikipedia.org/wiki/Multiple_dispatch)
# mechanism. When an algorithm calls a function such as `getobs!(obs, env)`, Julia will select from
# all functions with that name depending on `typeof(obs)` and `typeof(env)`.
struct Humanoid{S<:MJSim} <: AbstractMuJoCoEnvironment
    sim::S
end

LyceumMuJoCo.getsim(env::Humanoid) = env.sim #src (needs to be here for below example to work)

# `Humanoid` (and all subtypes of `AbstractEnvironment`) are designed to be used in a single
# threaded context. To use `Humanoid` in a multi-threaded context, one could simply create
# `Threads.nthreads()` instances of `Humanoid`:
modelpath = joinpath(@__DIR__, "humanoid.xml")
envs = [Humanoid(MJSim(modelpath, skip = 2)) for i = 1:Threads.nthreads()]
Threads.@threads for i = 1:Threads.nthreads()
    thread_env = envs[Threads.threadid()]
    step!(thread_env)
end

# As `Humanoid` only ever uses its internal `jlModel` (found at `sim.m`) in a read-only
# fashion, we can make a performance optimization by sharing a single instance of `jlModel`
# across each thread, resulting in improved cache efficiency. `LyceumMuJoCo.tconstruct`,
# short for "thread construct", helps us to do just that by providing a common interface
# for defining "thread-aware" constructors. Below, we make a call to
# `tconstruct(MJSim, n, modelpath, skip = 2)` which will construct `n` instances of `MJSim`
# constructed from `modelpath` and with a `skip` of 2, all sharing the exact same `jlModel`
# instance, and return `n` instances of `Humanoid`. All of the environments provided by
# LyceumMuJoCo feature similar definitions of `tconstruct` as found below.
Humanoid() = first(tconstruct(Humanoid, 1))
function LyceumMuJoCo.tconstruct(::Type{Humanoid}, n::Integer)
    modelpath = joinpath(@__DIR__, "humanoid.xml")
    return Tuple(Humanoid(s) for s in tconstruct(MJSim, n, modelpath, skip = 2))
end

# We can then use `tconstruct` as follows:
envs = tconstruct(Humanoid, Threads.nthreads())
Threads.@threads for i = 1:Threads.nthreads()
    thread_env = envs[Threads.threadid()]
    step!(thread_env)
end

#md # ### Utilities

# The following are helpers for the "get-up" task we'd like to consider.
# We want the humanoid to stand up, thus we need to grab the model's height, as
# well as record a laying down position that we can use to set the state to.
# By exploring the model in the REPL or MJCF/XML file we can see that `sim.d.qpos[3]`
# is the index for the z-axis (height) of the root joint. The `LAYING_QPOS` data was
# collected externally by posing the model into a supine pose; one can either use
# LyceumMuJoCoViz or simulate.cpp included with a MuJoCo release to do this.
_getheight(shapedstate::ShapedView, ::Humanoid) = shapedstate.qpos[3]
const LAYING_QPOS = [
    -0.164158,
    0.0265899,
    0.101116,
    0.684044,
    -0.160277,
    -0.70823,
    -0.0693176,
    -0.1321,
    0.0203937,
    0.298099,
    0.0873523,
    0.00634907,
    0.117343,
    -0.0320319,
    -0.619764,
    0.0204114,
    -0.157038,
    0.0512385,
    0.115817,
    -0.0320437,
    -0.617078,
    -0.00153819,
    0.13926,
    -1.01785,
    -1.57189,
    -0.0914509,
    0.708539,
    -1.57187,
];

#md # ### The `AbstractMuJoCoEnvironment` and `AbstractEnvironment` APIs

# LyceumMuJoCo requires access to the underlying `MJSim` simulator, thus any LyceumMuJoCo
# environments need to point to the correct field in the environment struct that is the
# simulator; in our case there's only one field: `sim`.
LyceumMuJoCo.getsim(env::Humanoid) = env.sim

# Normally we could rely on MuJoCo to reset the model to the default configuration when the
# model XML is loaded; the humanoid.xml model, however, defaults to a vertical pose. To
# reset the model to our laying down or supine pose, we can copy in the data from
# `LAYING_QPOS` above to `d.qpos`. Calling `forward!` here is the same as
# `mj_forward(env.sim.m, env.sim.d)`, for a pure MuJoCo reference.
function LyceumMuJoCo.reset!(env::Humanoid)
    reset!(env.sim)
    env.sim.d.qpos .= LAYING_QPOS
    forward!(env.sim)
    return env
end

# This reward function uses the `_getheight` helper function above to get the model's
# height when the function is called. We also specify a target height of 1.25 and
# penalize the agent for deviating from the target height. There is also a small penalty
# for using large control activations; if the coefficient is made larger, the agent
# may not move at all!
function LyceumMuJoCo.getreward(state, action, obs, env::Humanoid)
    height = _getheight(statespace(env)(state), env)
    target = 1.25
    reward = 1.0
    if height < target
        reward -= 2.0 * abs(target - height)
    end
    reward -= 1e-3 * norm(action)^2

    return reward
end

# Finally, we can specify an evaluation function. The difference between `geteval` and
# `getreward` is that `getreward` is the shaped reward our algorithm is optimizing for,
# while `geteval` lets us track a useful value for monitoring performance, such as height.
# Plotting this eval function will show the agent's height over time and is very useful
# for reviewing actual desired behavior, regardless of the reward achieved, as it can be
# used to diagnose reward specification problems.
function LyceumMuJoCo.geteval(state, action, obs, env::Humanoid)
    return _getheight(statespace(env)(state), env)
end

#md # ### Running Experiments

# As discussed in the [Julia performance tips](https://docs.julialang.org/en/v1/manual/performance-tips/),
# globals can hinder performance. To avoid this, we construct the `MPPI` and
# `ControllerIterator` instances within a function. This also lets us easily run our
# experiment with different choices of parameters (e.g. `H`). Like most algorithms in
# `LyceumAI`, `MPPI` accepts a "thread-aware" environment constructor as well as any
# algorithm parameters. In this case, we just pass a closure around the `tconstruct`
# function we defined above. `MPPI`, being a single-step algorithm, is itself not iterable,
# so we wrap it in a `ControllerIterator` which simply calls the passed-in function
# `ctrlfn` for `T` timesteps, while simultaneously plotting and logging the trajectory rollout.
function humanoid_MPPI(etype = Humanoid; T = 200, H = 64, K = 64)
    env = etype()

    ## The following parameters work well for this get-up task, and may work for
    ## similar tasks, but are not invariant to the model.
    mppi = MPPI(
        env_tconstructor = n -> tconstruct(etype, n),
        covar = Diagonal(0.05^2 * I, size(actionspace(env), 1)),
        lambda = 0.4,
        H = H,
        K = K,
        gamma = 1.0,
    )

    ctrlfn = (action, state, obs) -> getaction!(action, state, mppi)
    iter = ControllerIterator(ctrlfn, env; T = T, plotiter = div(T, 10))

    ## We can time the following loop; if it ends up less than the time the
    ## MuJoCo models integrated forward in, then one could conceivably run this
    ## MPPI MPC controller interactively...
    elapsed = @elapsed for (t, traj) in iter
        ## If desired, one can inspect `traj`, `env`, or `mppi` at each timestep.
    end

    if elapsed < time(env)
        @info "We ran in real time!"
    end

    ## Save our experiment results to a file for later review.
    savepath = "/tmp/opt_humanoid.jlso"
    exper = Experiment(savepath, overwrite = true)
    exper[:etype] = etype

    for (k, v) in pairs(iter.trajectory)
        exper[k] = v
    end
    finish!(exper)

    return mppi, env, iter.trajectory
end

#md # ### Checking Results

seed_threadrngs!(1) #src

# The MPPI algorithm, and any that you develop, can and should use plotting tools
# to track progress as they go.
mppi, env, traj = humanoid_MPPI();
plot(
    [traj.rewards traj.evaluations],
    labels = ["Reward" "Evaluation"],
    title = "Humanoid Standup",
    legend = :bottomright,
)

# If one wanted to review the results after training, or prepare plots for presentations,
# one can load the data from disk instead.
data = JLSO.load("/tmp/opt_humanoid.jlso")
plot(
    [data["rewards"] data["evaluations"]],
    labels = ["Reward" "Evaluation"],
    title = "Humanoid Standup",
    legend = :bottomright,
)

using Test #src
@test abs(geteval(env)) > 1.2 #src
@test data["rewards"] == traj.rewards #src


text:
include("Interval.jl")
include("PeriodicInterval.jl")
include("Ray.jl")
include("Circle.jl")
include("Line.jl")
include("Arc.jl")

include("UnionDomain.jl")
include("PiecewiseInterval.jl")
include("Curve.jl")

include("Point.jl")

typealias AffineDomain Union{Interval,PeriodicInterval,Ray,Line}

# These are needed for spaces to auto-convert [a,b] to Interval
function Base.convert{T<:Number}(::Type{Domain},d::AbstractVector{T})
    @assert length(d) >1

    if length(d) == 2
        if abs(d[1]) == Inf && abs(d[2]) == Inf
            Line(d)
        elseif abs(d[2]) == Inf || abs(d[1]) == Inf
            Ray(d)
        else
            Interval(d[1],d[2])
        end
    else
        # TODO: convert to PiecewiseInterval
        Domain(d[1:2]) ∪ Domain(d[2:end])
    end
end

function Base.convert{T<:Number}(::Type{PeriodicDomain},d::AbstractVector{T})
    @assert length(d) == 2

    if abs(d[1]) == Inf
        PeriodicLine(d)
    else
        PeriodicInterval(d[1],d[2])
    end
end

Base.promote_rule{D<:Domain,T<:Number}(::Type{D},::Type{Vector{T}})=UnivariateDomain{T}
Base.promote_rule{D<:PeriodicDomain,T<:Number}(::Type{D},::Type{Vector{T}})=PeriodicDomain{T}

#issubset between domains

Base.issubset(a::PeriodicInterval,b::Interval)=Interval(a.a,a.b)⊆b
Base.issubset(a::Interval,b::PeriodicInterval)=a⊆PeriodicInterval(b.a,b.b)
Base.issubset{T<:Real}(a::Interval{T},b::PiecewiseInterval{T})=a⊆Interval(first(b.points),last(b.points))
Base.issubset(a::Interval,b::Line)=first(a)∈b && last(a)∈b

function Base.intersect(a::Interval,b::Line)
    @assert a ⊆ b
    a
end

Base.intersect(b::Line,a::Interval)=intersect(a,b)

function Base.setdiff(b::Line,a::Interval)
    @assert a ⊆ b
    if first(a)>last(a)
        b\reverse(a)
    else
        Ray([first(b),first(a)]) ∪ Ray([last(a),last(b)])
    end
end

function Base.setdiff(b::Interval,a::Point)
    if !(a ⊆ b)
        a
    elseif first(b)==a.x  || last(b) == a.x
        a
    else
        Interval(first(b),a.x) ∪ Interval(a.x,last(b))
    end
end

# sort

Base.isless{T1<:Real,T2<:Real}(d1::Interval{T1},d2::Ray{false,T2})=d1≤d2.center
Base.isless{T1<:Real,T2<:Real}(d2::Ray{true,T2},d1::Interval{T1})=d2.center≤d1

# multivariate domainxs

include("multivariate.jl")
include("Disk.jl")


text:
using Test
using DiffEqPhysics, ForwardDiff, OrdinaryDiffEq
using StaticArrays, LinearAlgebra, Random

test_solve(prob...) = mapreduce(p->solve(p, VelocityVerlet(), dt=1//2).u, ==, prob)

H(dθ, θ, p) = dθ / 2 - 9.8 * cos(θ)
p0, q0   = rand(2)
acc      = (v, x, p, t) -> ForwardDiff.derivative(x -> -H(v[1], x, p), x[1])
vel      = (v, x, p, t) -> ForwardDiff.derivative(v -> H(v, x[1], p), v[1])
prob1    = HamiltonianProblem(H, p0, q0, (0., 10.))
prob_1   = DynamicalODEProblem(acc, vel, p0, q0, (0., 10.))
@test test_solve(prob1, prob_1)

p0 = @SVector rand(2)
q0 = @SVector rand(2)
H(dθ, θ, p) = dθ[1] / 2 + dθ[2] / 2 - 9.8 * cos(θ[1]) - 9.8 * cos(θ[2])
prob2    = HamiltonianProblem(H, p0, q0, (0., 10.))
acc      = (v, x, p, t) -> ForwardDiff.gradient(x -> -H(v, x, p), x)
vel      = (v, x, p, t) -> ForwardDiff.gradient(v -> H(v, x, p), v)
prob_2   = DynamicalODEProblem(acc, vel, p0, q0, (0., 10.))
@test test_solve(prob2, prob_2)

H(dθ, θ, p) = (dθ / 2 - 9.8 * cos.(θ))[1]
p0, q0 = [rand(5) for i in 1:2]
acc = (dv, v, x, p, t) -> ForwardDiff.gradient!(dv, x -> -H(v, x, p), x)
vel = (dx, v, x, p, t) -> ForwardDiff.gradient!(dx, v -> H(v, x, p), v)
prob3  = HamiltonianProblem(H, p0, q0, (0., 10.))
prob_3 = DynamicalODEProblem(acc, vel, p0, q0, (0., 10.))
@test test_solve(prob3, prob_3)


text:
#!/usr/bin/env julia
#
# Called with no arguments will create the system image
#     ClimateMachine.so
# in the `@__DIR__` directory.
#
# Called with a single argument the system image will be placed in the path
# specified by the argument (relative to the callers path)
sysimage_path =
    isempty(ARGS) ? joinpath(@__DIR__, "ClimateMachine.so") : ARGS[1]

using Pkg
Pkg.add("PackageCompiler")
using PackageCompiler
Pkg.activate(joinpath(@__DIR__, "..", ".."))
Pkg.instantiate()
pkgs = Pkg.installed()
delete!(pkgs, "Pkg")

# "TF_BUILD" is check to see if we are running on Azure for CI
if !haskey(ENV, "TF_BUILD")
    PackageCompiler.create_sysimage(
        [Symbol(s) for s in keys(pkgs)];
        sysimage_path = sysimage_path,
        precompile_execution_file = joinpath(
            @__DIR__,
            "..",
            "..",
            "test",
            "Numerics",
            "DGmethods",
            "Euler",
            "isentropicvortex.jl",
        ),
    )
else
    return true
end


text:
export PermMatrix, pmrand

"""
    PermMatrix{Tv, Ti}(perm::AbstractVector{Ti}, vals::AbstractVector{Tv}) where {Tv, Ti<:Integer}
    PermMatrix(perm::Vector{Ti}, vals::Vector{Tv}) where {Tv, Ti}
    PermMatrix(ds::AbstractMatrix)

PermMatrix represents a special kind linear operator: Permute and Multiply, which means `M * v = v[perm] * val`
Optimizations are used to make it much faster than `SparseMatrixCSC`.

* `perm` is the permutation order,
* `vals` is the multiplication factor.

[Generalized Permutation Matrix](https://en.wikipedia.org/wiki/Generalized_permutation_matrix)
"""
struct PermMatrix{Tv, Ti<:Integer, Vv<:AbstractVector{Tv}, Vi<:AbstractVector{Ti}} <: AbstractMatrix{Tv}
    perm::Vi   # new orders
    vals::Vv   # multiplied values.

    function PermMatrix{Tv, Ti, Vv, Vi}(perm::Vi, vals::Vv) where {Tv, Ti<:Integer, Vv<:AbstractVector{Tv}, Vi<:AbstractVector{Ti}}
        if length(perm) != length(vals)
            throw(DimensionMismatch("permutation ($(length(perm))) and multiply ($(length(vals))) length mismatch."))
        end
        new{Tv, Ti, Vv, Vi}(perm, vals)
    end
end

function PermMatrix{Tv, Ti}(perm, vals) where {Tv, Ti<:Integer}
    PermMatrix{Tv, Ti, Vector{Tv}, Vector{Ti}}(Vector{Ti}(perm), Vector{Tv}(vals))
end

function PermMatrix(perm::Vi, vals::Vv) where {Tv, Ti<:Integer, Vv<:AbstractVector{Tv}, Vi<:AbstractVector{Ti}}
    PermMatrix{Tv,Ti, Vv, Vi}(perm, vals)
end

################# Array Functions ##################

size(M::PermMatrix) = (length(M.perm), length(M.perm))
function size(A::PermMatrix, d::Integer)
    if d < 1
        throw(ArgumentError("dimension must be ≥ 1, got $d"))
    elseif d<=2
        return length(A.perm)
    else
        return 1
    end
end
getindex(M::PermMatrix, i::Integer, j::Integer) = M.perm[i] == j ? M.vals[i] : 0
copyto!(A::PermMatrix, B::PermMatrix) = (copyto!(A.perm, B.perm); copyto!(A.vals, B.vals); A)

"""
    pmrand(T::Type, n::Int) -> PermMatrix

Return random PermMatrix.
"""
function pmrand end

pmrand(::Type{T}, n::Int) where T = PermMatrix(randperm(n), randn(T, n))
pmrand(n::Int) = pmrand(Float64, n)

similar(x::PermMatrix{Tv, Ti}) where {Tv, Ti} = PermMatrix{Tv, Ti}(similar(x.perm), similar(x.vals))
similar(x::PermMatrix{Tv, Ti}, ::Type{T}) where {Tv, Ti, T} = PermMatrix{T, Ti}(similar(x.perm), similar(x.vals, T))

# TODO: rewrite this
# function show(io::IO, M::PermMatrix)
#     println("PermMatrix")
#     for item in zip(M.perm, M.vals)
#         i, p = item
#         println("- ($i) * $p")
#     end
# end

######### sparse array interfaces  #########
nnz(M::PermMatrix) = length(M.vals)
nonzeros(M::PermMatrix) = M.vals
dropzeros!(M::PermMatrix; trim::Bool=false) = M
isdense(::PermMatrix) = false

Base.zero(pm::PermMatrix) = PermMatrix(pm.perm, zero(pm.vals))


text:
function ueberpruefePixelkanten(X,Y)
if X == [-1.5 -0.5 0.5 1.5]
    if Y == [-1.5 -0.5 0.5 1.5]
        println("Ihr habt die Pixelkanten korrekt gewählt!")
    else
        println("Schaut euch noch einmal die Pixelkanten bei Y an.")
    end
else
    println("Überprüft noch einmal eure Eingaben.")
end
end


text:
using MathOptInterface
const MOI = MathOptInterface

# Some tests are excluded because UniversalFallback accepts absolutely
# everything.

MOI.Test.runtests(
    MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    ),
    MOI.Test.Config(),
    exclude = [
        "test_model_ScalarFunctionConstantNotZero",
        "test_model_copy_to_UnsupportedAttribute",
        "test_model_copy_to_UnsupportedConstraint",
        "test_model_supports_constraint_ScalarAffineFunction_EqualTo",
        "test_model_supports_constraint_VariableIndex_EqualTo",
        "test_model_supports_constraint_VectorOfVariables_Nonnegatives",
    ],
    warn_unsupported = true,
)

# Run the previously excluded tests, this time without UniversalFallback.

MOI.Test.runtests(
    MOI.Utilities.MockOptimizer(
        MOI.Utilities.Model{Float64}(),
        scalar_function_constant_non_zero = true,
    ),
    MOI.Test.Config(),
    include = [
        "test_model_ScalarFunctionConstantNotZero",
        "test_model_copy_to_UnsupportedAttribute",
        "test_model_copy_to_UnsupportedConstraint",
        "test_model_supports_constraint_ScalarAffineFunction_EqualTo",
        "test_model_supports_constraint_VariableIndex_EqualTo",
        "test_model_supports_constraint_VectorOfVariables_Nonnegatives",
    ],
)


text:
# open up a block to limit scope and lifetime of everything to
# ensure no cross-talk between tests
begin

    N = 500

    A = rand(N)
    B = rand(N)
    C = rand(N)

    # addition
    @et add_et = A + B + C
    add_native = A + B + C

    @test add_et == add_native;

    #subtraction
    @et sub_et = A - B - C
    sub_native = A - B - C

    @test sub_et == sub_native;

    # elementwise mult
    @et dottimes_et = A .* B .* C
    dottimes_native = A .* B .* C

    @test dottimes_et == dottimes_native;

    # elementwise div
    @et dotdiv_et = A ./ B ./ C
    dotdiv_native = A ./ B ./ C

    @test dotdiv_et == dotdiv_native;

end

begin

    N = 50

    A = rand(N,N)
    B = rand(N,N)
    C = rand(N,N)

    #test order-2 arrays
    @et res_et = A+B+C
    res_native = A+B+C

    @test res_et == res_native

end

begin

    N = 15

    A = rand(N,N,N)
    B = rand(N,N,N)
    C = rand(N,N,N)

    #test order-3 arrays
    @et res_et = A+B+C
    res_native = A+B+C

    @test res_et == res_native

    # test .* and ./
    @et res_et = A.*B./C
    res_native = A.*B./C

    @test res_et == res_native
end


text:
mutable struct MultivariateGaussianDist
    _mu::AbstractArray{Float64,1}
    _cov::AbstractArray{Float64,2}
    function MultivariateGaussianDist(dim::Int64)
        new(zeros(dim), Matrix{Float64}(I, dim, dim))
    end
end

function pdf(gaussian::MultivariateGaussianDist, x::AbstractArray{Float64,1})
    @assert size(x)[1] == size(gaussian._mu)[1]
    N = size(x)[1]
    d = x - gaussian._mu
    quadratics = transpose(d) * gaussian._cov * d
    @assert typeof(quadratics) == Float64
    return exp(-0.5 * quadratics) / sqrt(det(gaussian._cov)) / (2 * pi)^(0.5 * N)
end

function pdf(gaussian::MultivariateGaussianDist, X::AbstractArray{Float64,2})
    # X consists of list of samples [x_1, x_2, ,,, x_(n_samples)], and x_i is a vertical vector
    @assert size(X)[1] == size(gaussian._mu)[1]
    N = size(X)[1]
    n_samples = size(X)[2]
    # In Julia lang, 1-D AbstractArray is vertical vector
    d = X .- gaussian._mu
    #quadratics = diag(transpose(d) * gaussian._cov * d);
    quadratics = [transpose(d[:, i]) * gaussian._cov * d[:, i] for i = 1:n_samples]
    @assert ndims(quadratics) == 1 && size(quadratics)[1] == n_samples
    return exp.(-0.5 * quadratics) / sqrt(det(gaussian._cov)) / (2 * pi)^(0.5 * N)
end

function fitting(gaussian::MultivariateGaussianDist, X::AbstractArray{Float64,1})
    gaussian._mu = X
    gaussian._cov = cov(X, corrected = true)
end

function fitting(gaussian::MultivariateGaussianDist, X::AbstractArray{Float64,2})
    # X consists of the list of samples [x_1, x_2, ,,, x_(n_samples)]
    @assert size(X)[1] == size(gaussian._mu)[1]
    n_samples = size(X)[2]
    N = size(X)[1]
    gaussian._mu = collect(reshape(mean(X, dims = 2), N))
    gaussian._cov = cov(X, dims = 2, corrected = true)
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
# Bounding box
# ---------------------------------------------------------------------------------------------
# Bounding box of a vector of points
function boundingbox(points::Vector{<:Point2D})
    xmin = ymin = typemax(T)
    xmax = ymax = typemin(T)
    for i = 1:length(points)
        x,y = points[i].coord
        if x < xmin
            xmin = x
        end
        if xmax < x
            xmax = x
        end
        if y < ymin
            ymin = y
        end
        if ymax < y
            ymax = y
        end
    end
    return AABox2D(Point2D(xmin, ymin),
                   Point2D(xmax, ymax))
end

function boundingbox(points::SVector{L, Point2D{T}}) where {L,T}
    xmin = ymin = typemax(T)
    xmax = ymax = typemin(T)
    for i = 1:L
        x,y = points[i].coord
        if x < xmin
            xmin = x
        end
        if xmax < x
            xmax = x
        end
        if y < ymin
            ymin = y
        end
        if ymax < y
            ymax = y
        end
    end
    return AABox2D(Point2D(xmin, ymin),
                   Point2D(xmax, ymax))
end

# Bounding box of a vector of points
function boundingbox(points::Vector{<:Point3D})
    xmin = ymin = zmin = typemax(T)
    xmax = ymax = zmax = typemin(T)
    for i = 1:length(points)
        x,y,z = points[i].coord
        if x < xmin
            xmin = x
        end
        if xmax < x
            xmax = x
        end
        if y < ymin
            ymin = y
        end
        if ymax < y
            ymax = y
        end
        if z < zmin
            zmin = z
        end
        if zmax < z
            zmax = z
        end
    end
    return AABox3D(Point3D(xmin, ymin, zmin),
                   Point3D(xmax, ymax, zmax))
end

function boundingbox(points::SVector{L, Point3D{T}}) where {L,T}
    xmin = ymin = zmin = typemax(T)
    xmax = ymax = zmax = typemin(T)
    for i = 1:L
        x,y,z = points[i].coord
        if x < xmin
            xmin = x
        end
        if xmax < x
            xmax = x
        end
        if y < ymin
            ymin = y
        end
        if ymax < y
            ymax = y
        end
        if z < zmin
            zmin = z
        end
        if zmax < z
            zmax = z
        end
    end
    return AABox3D(Point3D(xmin, ymin, zmin),
                   Point3D(xmax, ymax, zmax))
end

# Axis-aligned bounding box
# ---------------------------------------------------------------------------------------------
# Find the axis-aligned bounding box of the segment
#
# Find the extrema for x and y by finding the r_x such that dx/dr = 0
# and r_y such that dy/dr = 0
# 𝗾(r) = r²𝘂 + r𝘃 + 𝘅₁
# 𝗾′(r) = 2r𝘂 + 𝘃 ⟹  r_x, r_y = -𝘃 ./ 2𝘂
# Compare the extrema with the segment's endpoints to find the AABox
function boundingbox(q::QuadraticSegment{N}) where {N}
    𝘂 = q.𝘂
    𝘃 = q.𝘃
    𝗿 = 𝘃 ./ -2𝘂
    𝗽_stationary = 𝗿*𝗿*𝘂 + 𝗿*𝘃 + q.𝘅₁
    𝗽_min = min.(q.𝘅₁.coord, q.𝘅₂.coord)
    𝗽_max = max.(q.𝘅₁.coord, q.𝘅₂.coord)
    if N === 2
        xmin, ymin = 𝗽_min
        xmax, ymax = 𝗽_max
        if 0 < 𝗿[1] < 1
            xmin = min(𝗽_min[1], 𝗽_stationary[1])
            xmax = max(𝗽_max[1], 𝗽_stationary[1])
        end
        if 0 < 𝗿[2] < 1
            ymin = min(𝗽_min[2], 𝗽_stationary[2])
            ymax = max(𝗽_max[2], 𝗽_stationary[2])
        end
        return AABox2D(Point2D(xmin, ymin), Point2D(xmax, ymax))
    else # N === 3
        xmin, ymin, zmin = 𝗽_min
        xmax, ymax, zmax = 𝗽_max
        if 0 < 𝗿[1] < 1
            xmin = min(𝗽_min[1], 𝗽_stationary[1])
            xmax = max(𝗽_max[1], 𝗽_stationary[1])
        end
        if 0 < 𝗿[2] < 1
            ymin = min(𝗽_min[2], 𝗽_stationary[2])
            ymax = max(𝗽_max[2], 𝗽_stationary[2])
        end
        if 0 < 𝗿[3] < 1
            zmin = min(𝗽_min[3], 𝗽_stationary[3])
            zmax = max(𝗽_max[3], 𝗽_stationary[3])
        end
        return AABox3D(Point3D(xmin, ymin, zmin), Point3D(xmax, ymax, zmax))
    end
end

# Return the AABox which contains both bb₁ and bb₂
function Base.union(bb₁::AABox{Dim, T}, bb₂::AABox{Dim, T}) where {Dim, T}
    return AABox(Point{Dim, T}(min.(bb₁.minima.coord, bb₂.minima.coord)),
                 Point{Dim, T}(max.(bb₁.maxima.coord, bb₂.maxima.coord)))
end

# Bounding box
# ---------------------------------------------------------------------------------------------
boundingbox(poly::Polygon) = boundingbox(poly.points)

# Axis-aligned bounding box
function boundingbox(mesh::LinearUnstructuredMesh)
    # The bounding box may be determined entirely from the points.
    return boundingbox(mesh.points)
end

# Axis-aligned bounding box
function boundingbox(mesh::QuadraticUnstructuredMesh)
    return mapreduce(x->boundingbox(x), union, materialize_edges(mesh))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
######### CmdStan program example  ###########

using CmdStan

ProjDir = @__DIR__
cd(ProjDir) do

  bernoullimodel = "
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
    y ~ bernoulli(theta);
  }
  ";

  observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

  # Preserve these variables outside cd/do block.
  global stanmodel, rc, samples, cnames, df

  stanmodel = Stanmodel(name="bernoulli", model=bernoullimodel,
    printsummary=false);

  rc, samples, cnames = stan(stanmodel, observeddata, ProjDir, CmdStanDir=CMDSTAN_HOME);

  if rc == 0
    # Fetch cmdstan summary as a DataFrame
    df = read_summary(stanmodel)

    df |> display
    println()

    df[df.parameters .== :theta, [:mean, :ess, :r_hat]] |> display
  end

end # cd block


text:
# canonical data version 4.0.0

using Test

include("high-scores.jl")

@testset "List of scores" begin
    @test scores([30, 50, 20, 70]) == [30, 50, 20, 70]
end

@testset "Latest score" begin
    @test scores([100, 0, 90, 30], latest=true) == 30
end

@testset "Personal best" begin
    @test scores([40, 100, 70], top=1) == 100
end

@testset "Top 3 scores" begin
    @test scores([10, 30, 90, 30, 100, 20, 10, 0, 30, 40, 40, 70, 70], top=3) == [100, 90, 70]
end

@testset "Personal top highest to lowest" begin
    @test scores([20, 10, 30], top=3) == [30, 20, 10]
end

@testset "Personal top when there is a tie" begin
    @test scores([40, 20, 40, 30], top = 3) == [40, 40, 30]
end

@testset "Personal top when there are less than 3" begin
    @test scores([30, 70], top=3) == [70, 30]
end

@testset "Personal top when there is only one" begin
    @test scores([40], top=3) == [40]
end

max_stars_repo_name: josobar/julia

max_issues_repo_name: josobar/julia

max_forks_repo_name: josobar/julia