text:
using DifferentialEquations
using NLsolve
using Plots

# The range of values α = M₂/(M₁+M₂) where the Lagrange points will be found
αvals = 0.001:0.001:0.5
# The number of values of α
nαvals = length(αvals)
# Will store the positions of the Lagrange points
lagrangepoints = zeros(Float64, nαvals, 10);

for i ∈ 1:nαvals
    # Set the value of α = M₂/(M₁+M₂)
	α = αvals[i]

    # Set up the function for NLsolve to use
	function evalF!(F, X)
		x, y = X
		F[1] = x -
			(1 - α)*(x + α)/√((x + α)^2 + y^2)^3 -
			α*(x + α - 1)/√((x + α - 1)^2 + y^2)^3
		F[2] = y -
			(1 - α)*y/√((x + α)^2 + y^2)^3 -
			α*y/√((x + α - 1)^2 + y^2)^3
	end

    # Run nlsolve at each of the five predetermined starting positions to find the five Lagrange points
	lagrangepoints[i,1:2 ] = nlsolve(evalF!, [ 0.0, 0.0], autodiff=:forward).zero
	lagrangepoints[i,3:4 ] = nlsolve(evalF!, [ 1.0, 0.0], autodiff=:forward).zero
	lagrangepoints[i,5:6 ] = nlsolve(evalF!, [-1.0, 0.0], autodiff=:forward).zero
	lagrangepoints[i,7:8 ] = nlsolve(evalF!, [ 0.0, 1.0], autodiff=:forward).zero
	lagrangepoints[i,9:10] = nlsolve(evalF!, [ 0.0,-1.0], autodiff=:forward).zero

end

# Plot the values of α against the coordinates of the Lagrange points,
#  both as found numerically, and as approximated theoretically
function locationplots()
    return plot(
        plot(αvals, [lagrangepoints[:,1 ], α -> 1 - ∛(α/3)]),
        plot(αvals, [lagrangepoints[:,2 ], α -> 0]),
        plot(αvals, [lagrangepoints[:,3 ], α -> 1 + ∛(α/3)]),
        plot(αvals, [lagrangepoints[:,4 ], α -> 0]),
        plot(αvals, [lagrangepoints[:,5 ], α -> - 1 - 5α/12]),
        plot(αvals, [lagrangepoints[:,6 ], α -> 0]),
        plot(αvals, [lagrangepoints[:,7 ], α -> 1/2 - α]),
        plot(αvals, [lagrangepoints[:,8 ], α -> (√3)/2]),
        plot(αvals, [lagrangepoints[:,9 ], α -> 1/2 - α]),
        plot(αvals, [lagrangepoints[:,10], α -> -(√3)/2]),
        linestyle = [:solid :dash],
        layout = (5,2),
        size = (600, 1000),
        label = ["numerical" "theoretical"],
        yticks = false,
        title = ["L₁ x" "L₁ y" "L₂ x" "L₂ y" "L₃ x" "L₃ y" "L₄ x" "L₄ y" "L₅ x" "L₅ y"]
    )
end

# Calculates the acceleration of m, given position r and velocity v, with p a vector of parameters
function acceleration(r, v, p)
    # Parameters for the system
    G, M₁, M₂, r₁, r₂, Ω = p

    # The four forces (including fictitious) acting on m
    gravitational₁ = - ( G * M₁ / hypot((r - r₁)...)^3 ) .* (r - r₁)
    gravitational₂ = - ( G * M₂ / hypot((r - r₂)...)^3 ) .* (r - r₂)
    coriolis = - 2Ω .* [[0, -1] [1, 0]] * v
    centrifugal = (Ω^2) .* r
    return gravitational₁ + gravitational₂ + coriolis + centrifugal
end

# Sets up the derivative function for DifferentialEquations to use
function derivative!(du, u, p, t)
    du[1] = u[3]
    du[2] = u[4]
    du[3], du[4] = acceleration(u[1:2], u[3:4], p)
end

# Uses DifferentialEquations to simulate the trajectory of m for time T with initial conditions
#  r₀ and v₀, and parameters α and M₁
function simulate(α, G, r₀, v₀, T)
    # M₂ is the mass of the other large stationary body
    M₂ = M₁ * α/(1-α)

    # r₁ and r₂ are the positions of M₁ and M₂
    r₁ = [-α, 0]
    r₂ = [1-α, 0]

	# G is the gravitational constant, which is 1 for an appropriate system of units
	G = 1
    # Ω is the magnitude of the angular velocity of the rotating frame of reference
    Ω = √(G * (M₁ + M₂))

    # DifferentialEquations simulates the mass for time T
    return solve(ODEProblem(derivative!, vcat(r₀, v₀), (0.0, float(T)), [G, M₁, M₂, r₁, r₂, Ω]))
end

# Creates a diagram of the Sun/Earth/Moon system to test the simulation
function moondiagram()
    diagram = plot(
		simulate(3e-6, 39, [1.003, 0.0], [0.0, 0.215], 1/12),
		vars = (1,2),
		linecolor = :white,
		linewidth = 1,
		arrow = true,
		bg = :black,
		ticks = false,
		showaxis = false,
		legend = false,
		xlims = [0.995, 1.005],
		ylims = [-0.005, 0.005]
	)

	scatter!(
		diagram,
		[(1-(3e-6),0)],
		markercolor = :aqua,
		markersize = 4
	)

    return diagram
end

# Plots the trajectory of m starting near Lₙ up to time T, with parameters α and M₁
function trajectory(α, M₁, n, T, lims = [2,1])
	i = findfirst(≥(α), αvals)
    simulatedtrajectory = simulate(αvals[i], M₁, lagrangepoints[i, (2n-1):2n], [0.0, 0.0], T)

	diagram = plot(
		simulatedtrajectory,
		vars = (1,2),
		linecolor = :white,
		linewidth = 1,
		bg = :black,
		ticks = false,
		showaxis = false,
		legend = false,
        xlims = [-lims[1],lims[1]],
        ylims = [-lims[2],lims[2]],
        size = Tuple(lims) .* round(Int64, 1200/sum(lims))
	)

	scatter!(
		diagram,
		[Tuple(simulatedtrajectory.u[1][1:2]), (-α,0), (1-α, 0)],
		markercolor = [:white, :yellow, :aqua],
		markersize = [3, 10, 4]
	)

	return diagram
end


text:
using PyPlot
using PyCall
@pyimport matplotlib.animation as anim

srand(12345)

b = 0.8 # 相関項(Correlation term)
x = 10.0
y = 3.0
times = 200
ims = []

for i = 1:times
    u1 = rand()
    u2 = rand()
    x_old = x
    x = sqrt(-2*log(u1))*cos(2*pi*u2) + b*y
    push!(ims, ([x_old, x], [y, y])) # 表示せずに配列imsに保存
    y_old = y
    y = sqrt(-2*log(u1))*sin(2*pi*u2) + b*x
    push!(ims, ([x, x], [y_old, y])) # 表示せずに配列imsに保存
end

########## display file
function displayfile(mimetype, filename)
    open(filename) do f
        base64text = base64encode(f)
        display("text/html", """<img src="data:$mimetype;base64,$base64text">""")
    end
end

########## initialize the frame
xs = [ims[i][j][1] for j in 1:2, i in 1:length(ims)]
ys = [ims[i][j][2] for j in 1:2, i in 1:length(ims)]
xmin, xmax = minimum(xs)-0.5, maximum(xs)+0.5
ymin, ymax = minimum(ys)-0.5, maximum(ys)+0.5
function initframe()
    xlim(xmin, xmax)
    ylim(ymin, ymax)
    grid(ls=":")
end

########## update the frame
function updateframe(t)
    # clf() #←全部消して書き直す場合にはコメントを外す
    #---------- フレームのアップデート開始
    if t ≥ 1
        X, Y = ims[t]
        plot(X, Y, lw=1.0)
    end
    title("t = $t")
    #---------- フレームのアップデート終了
    # plot() #←念のためのおまじない. []と書いてもよい. いらない場合も多い.
end

########## filename
filename = "gibbs_sampling.gif"

######### Construct Figure and Plot Data
fig = figure(figsize=(6.4, 4.8))

# Create the animation object by anim.FuncAnimaton
frames = collect(0:length(ims))
@time myanim = anim.FuncAnimation(fig, updateframe, init_func=initframe, frames=frames, interval=100)

# Convert it to an animated GIF file by Imagemagick.
@time myanim[:save](filename, writer="imagemagick")

sleep(0.1)

# Display the movie.
displayfile("image/gif", filename)

# Don't display the axes of figure().
clf()

ext: jl

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
name = gensym()
@eval D @discord_object struct $name
    a::Snowflake
    b::Int
    c::String
end

T = getfield(D, name)
@test fieldnames(T) == (:a, :b, :c)
@test map(n -> fieldtype(T, n), fieldnames(T)) == (
    Union{D.Snowflake, D.A.Null, Nothing, Missing},
    Union{Int, D.A.Null, Nothing, Missing},
    Union{String, D.A.Null, Nothing, Missing},
)

x = T()
@test x.a === missing && x.b === missing && x.c === missing
x = JSON3.read("{}", T)
@test x.a === missing && x.b === missing && x.c === missing
@test JSON3.write(x) == "{}"
x = JSON3.read("""{"b":10}""", T)
@test x.a === missing && x.b == 10 && x.c === missing
@test JSON3.write(x) == """{"b":10}"""
x = JSON3.read("""{"b":null}""", T)
@test x.a === missing && x.b === nothing && x.c === missing
@test JSON3.write(x) == """{"b":null}"""


text:
#   Unit examples of the FourierAnalysis Package for julia language
#
#   MIT License
#   Copyright (c) 2019-2022,
#   Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This example shows how to compute univariate measures in
#   the time-frequency domain.

#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #
#                                                                             #
#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #

using FourierAnalysis, FFTW, LinearAlgebra, Statistics, Plots

# add module for reading the two EEG text files to be used ater
push!(LOAD_PATH, @__DIR__)
using IOtxt

# Get EEG file names with complete path (they have extension .txt)
S=getFilesInDir(@__DIR__; ext=(".txt",))

# read the two EEG data files and put them in a Matrix object
X1=readEEG(S[1])
X2=readEEG(S[2])

#############################################################################
# Bivariate measures
# see Congedo, 2018: https://hal.archives-ouvertes.fr/hal-01868538/document

Pz=15
Fz=5
𝐱₁=[X1[:, Pz], X2[:, Pz]] # get the two times-series at electrode Pz
𝐱₂=[X1[:, Fz], X2[:, Fz]] # get the two times-series at electrode Fz
sr, wl, bandwidth=128, 512, 2
𝐘₁=TFanalyticsignal(𝐱₁, sr, wl, bandwidth; fmax=32, nonlinear=false)
𝐘₂=TFanalyticsignal(𝐱₂, sr, wl, bandwidth; fmax=32, nonlinear=false)
𝐀₁=TFamplitude(𝐘₁)
𝐀₂=TFamplitude(𝐘₂)

# Comodulation Eq. 0.11 (Com)
# compute the Com averaging in a TF region from a TFAnalyticSignalVector object
# NB If you compute it from Analytic Signal or Amplitude objects, those
# must be linear (note the nonlinear=false above)
Com=comodulation(𝐘₁, 𝐘₂, (8, 12), (1, 512); mode=mean)
# compute the Com averaging in a TF region from a TFAmplitudeVector object
Com=comodulation(𝐀₁, 𝐀₂, (8, 12), (1, 512); mode=mean)
# compute the Com averaging in a TF region directly from data
# In this care you don't have to worry about linearity of the analytic signal
Com=comodulation(𝐱₁, 𝐱₂, sr, wl, (8, 12), (1, 512), bandwidth; mode=mean)
# You can compute comodulation from smoothed amplitude:
Com=comodulation(𝐱₁, 𝐱₂, sr, wl, (8, 12), (1, 512), bandwidth;
                 mode=mean,
                 fsmoothing=blackmanSmoother,
                 tsmoothing=blackmanSmoother)

# you can go faster pre-computing a FFTW plan.
# This is aslo useful when you have to call the comodulation function several times
plan=Planner(plan_patient, 5, wl, Float64, true)
Com=comodulation(𝐱₁, 𝐱₂, sr, wl, (8, 12), (1, 512), bandwidth; mode=mean, planner=plan)

# compute the Com in a TF region from a TFAnalyticSignalVector object
Com=comodulation(𝐘₁, 𝐘₂, (8, 12), (1, 512); mode=extract)
# compute the Com in a TF region from a TFAmplitudeVector object
Com=comodulation(𝐀₁, 𝐀₂, (8, 12), (1, 512); mode=extract)
# compute the Com in a TF region directly from data
Com=comodulation(𝐱₁, 𝐱₂, sr, wl, (8, 12), (1, 512), bandwidth; mode=extract)

# All these operations can be done also for coherence measures, for example
# avoid conflict with package DSP, which also exports 'coherence'
Coh=FourierAnalysis.coherence(𝐘₁, 𝐘₂, (8, 12), (1, 512); mode=mean)
Coh=FourierAnalysis.coherence(𝐘₁, 𝐘₂, (8, 12), (1, 512); mode=extract)

# Compute all 5 coherence types
Coh=FourierAnalysis.coherence(𝐘₁, 𝐘₂, (8, 12), (1, 512); mode=extract, allkinds=true)

# phase coherence (phase-locking value)
𝐘₁=TFanalyticsignal(𝐱₁, sr, wl, bandwidth; fmax=32, nonlinear=true)
𝐘₂=TFanalyticsignal(𝐱₂, sr, wl, bandwidth; fmax=32, nonlinear=true)
Coh=FourierAnalysis.coherence(𝐘₁, 𝐘₂, (8, 12), (1, 512); mode=mean, nonlinear=true)

# or directly from data (no need to compute non-linear analytic signal in this case)
Coh=FourierAnalysis.coherence(𝐱₁, 𝐱₂, sr, wl, (8, 12), (1, 512), bandwidth; mode=mean, nonlinear=true)

# and also for non-linear meausures
# compute non-linear analyticSignal
𝐘₁=TFanalyticsignal(𝐱₁, sr, wl, bandwidth; fmax=32, nonlinear=true)
𝐘₂=TFanalyticsignal(𝐱₂, sr, wl, bandwidth; fmax=32, nonlinear=true)
# although you are allowed to compute the amplitude of non-linear
# analytic sygnal this way, it does not make sense as the amplitude
# is 1.0 everywhere, as you can check
𝐀₁=TFamplitude(𝐘₁)
𝐀₂=TFamplitude(𝐘₂)

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


text:
"""
### Nuclear attraction term
"""
function Aterm(i::Int64,r::Int64,u::Int64,l1::Int64,l2::Int64,ax::Float64,bx::Float64,cx::Float64,gamma::Float64)
  term1 = (-1)^(i+u)*binomial_prefactor(i,l1,l2,ax,bx)
  term2 = factorial(i)*cx^(i-2r-2u)
  term3 = (1/4/gamma)^(r+u)/factorial(r)/factorial(u)/factorial(i-2r-2u)
  return term1*term2*term3
end

function Aarray(l1::Int64,l2::Int64,a::Float64,b::Float64,c::Float64,g::Float64)
  Imax = l1+l2+1
  A = zeros(Float64,Imax)
  for i in 0:(Imax-1)
    for r in 0:div(i,2)
      for u in 0:div(i-2r,2)
        I = i-2r-u+1
        A[I] += Aterm(i,r,u,l1,l2,a,b,c,g)
      end
    end
  end
  return A
end

function nuclear_attraction(aexpn::Float64,ax::Float64,ay::Float64,az::Float64,
              aI::Int64,aJ::Int64,aK::Int64,
              bexpn::Float64,bx::Float64,by::Float64,bz::Float64,
              bI::Int64,bJ::Int64,bK::Int64,
              cx::Float64,cy::Float64,cz::Float64)
  #print("na($aexpn,$ax,$ay,$az,$aI,$aJ,$aK,$bexpn,$bx,$by,$bz,$bI,$bJ,$bK,$cx,$cy,$cz)=")
  px,py,pz = gaussian_product_center(aexpn,ax,ay,az,bexpn,bx,by,bz)
  gamma = aexpn+bexpn
  rab2 = dist2(ax-bx,ay-by,az-bz)
  rcp2 = dist2(cx-px,cy-py,cz-pz)
  Ax = Aarray(aI,bI,px-ax,px-bx,px-cx,gamma)
  Ay = Aarray(aJ,bJ,py-ay,py-by,py-cy,gamma)
  Az = Aarray(aK,bK,pz-az,pz-bz,pz-cz,gamma)
  total = 0
  for I in 0:(aI+bI)
    for J in 0:(aJ+bJ)
      for K in 0:(aK+bK)
        total += Ax[I+1]*Ay[J+1]*Az[K+1]*Fgamma(I+J+K,rcp2*gamma)
      end
    end
  end
  val=-2pi*exp(-aexpn*bexpn*rab2/gamma)*total/gamma
  #println(val)
  #println((Ax,Ay,Az,rcp2*gamma,Fgamma(0,rcp2*gamma)))
  return val
end

function nuclear_attraction(a::PGBF,b::PGBF,cx::Float64,cy::Float64,cz::Float64)
  return a.norm*b.norm*nuclear_attraction(a.expn,a.x,a.y,a.z,a.I,a.J,a.K,
                      b.expn,b.x,b.y,b.z,b.I,b.J,b.K,cx,cy,cz)
end
nuclear_attraction(a::PGBF,b::PGBF,c::Atom) = c.atno*nuclear_attraction(a,b,c.x,c.y,c.z)
nuclear_attraction(a::PGBF,b::PGBF,m::Molecule) = sum([nuclear_attraction(a,b,c) for c in m.atomlist])

function Fgamma(m::Int64,x::Float64,SMALL::Float64=1e-12)
  #println("Fgamma($m,$x)")
  x = max(x,SMALL) # Evidently needs underflow protection
  return 0.5*x^(-m-0.5)*gammainc(m+0.5,x)
end

function gammainc(a::Float64,x::Float64)
  # This is the series version of gamma from pyquante. For reasons I don't get, it
  # doesn't work around a=1. This works alright, but is only a stopgap solution
  # until Julia gets an incomplete gamma function programmed
  if abs(a-1) < 1e-3
    println("Warning: gammainc_series is known to have problems for a ~ 1")
  end
  if x < (a+1.0)
    #Use the series representation
    gam,gln = gser(a,x)
  else
    #Use continued fractions
    gamc,gln = gcf(a,x)
    gam = 1-gamc
  end
  return exp(gln)*gam
end

function gser(a::Float64,x::Float64,ITMAX::Int64=100,EPS::Float64=3e-9)
  # Series representation of Gamma. NumRec sect 6.1.
  gln=lgamma(a)
  if x == 0
    return 0,gln
  end
  ap = a
  delt = s = 1/a
  for i in 1:ITMAX
    ap += 1
    delt *= (x/ap)
    s += delt
    if abs(delt) < abs(s)*EPS
      break
    end
  end
  return s*exp(-x+a*log(x)-gln),gln
end

function gcf(a::Float64,x::Float64,ITMAX::Int64=200,EPS::Float64=3e-9,FPMIN::Float64=1e-30)
  #Continued fraction representation of Gamma. NumRec sect 6.1"
  gln=lgamma(a)
  b=x+1.-a
  c=1./FPMIN
  d=1./b
  h=d
  for i in 1:ITMAX
    an=-i*(i-a)
    b=b+2.
    d=an*d+b
    if abs(d) < FPMIN
      d=FPMIN
    end
    c=b+an/c
    if abs(c) < FPMIN
      c=FPMIN
    end
    d=1./d
    delt=d*c
    h=h*delt
    if abs(delt-1.) < EPS
      break
    end
  end
  gammcf = exp(-x+a*log(x)-gln)*h
  return gammcf,gln
end

# Need a nested scope to squeeze this into the contract function
function nuclear_attraction(a::CGBF,b::CGBF,cx::Float64,cy::Float64,cz::Float64)
  na(a,b) = nuclear_attraction(a,b,cx,cy,cz)
  contract(na,a,b)
end
function nuclear_attraction(a::CGBF,b::CGBF,c::Atom)
  na(a,b) = nuclear_attraction(a,b,c)
  contract(na,a,b)
end
function nuclear_attraction(a::CGBF,b::CGBF,m::Molecule)
  na(a,b) = nuclear_attraction(a,b,m)
  contract(na,a,b)
end

function test_a_terms()
  @assert Aterm(0,0,0,0,0,0.,0.,0.,0.) == 1.0
  @assert Aarray(0,0,0.,0.,0.,1.) == [1.0]
  @assert Aarray(0,1,1.,1.,1.,1.) == [1.0, -1.0]
  @assert Aarray(1,1,1.,1.,1.,1.) == [1.5, -2.5, 1.0]
  @assert Aterm(0,0,0,0,0,0.,0.,0.,1.) == 1.0
  @assert Aterm(0,0,0,0,1,1.,1.,1.,1.) == 1.0
  @assert Aterm(1,0,0,0,1,1.,1.,1.,1.) == -1.0
  @assert Aterm(0,0,0,1,1,1.,1.,1.,1.) == 1.0
  @assert Aterm(1,0,0,1,1,1.,1.,1.,1.) == -2.0
  @assert Aterm(2,0,0,1,1,1.,1.,1.,1.) == 1.0
  @assert Aterm(2,0,1,1,1,1.,1.,1.,1.) == -0.5
  @assert Aterm(2,1,0,1,1,1.,1.,1.,1.) == 0.5
end

function test_gamma()
  # gammainc test functions. Test values taken from Mathematica
  # println("a=0.5 test")
  @assert maximum([gammainc(0.5,float(x)) for x in 0:10]
      -Any[0, 1.49365, 1.69181, 1.7471, 1.76416, 1.76968,
        1.77151, 1.77213, 1.77234, 1.77241, 1.77244]) < 1e-5

  # println("a=1.5 test")
  @assert maximum([gammainc(1.5,float(x)) for x in 0:10]
      -Any[0, 1.49365, 1.69181, 1.7471, 1.76416, 1.76968,
        1.77151, 1.77213, 1.77234, 1.77241, 1.77244]) < 1e-5
  # println("a=2.5 test")
  @assert maximum([gammainc(2.5,float(x)) for x in 0:10]
      -Any[0, 0.200538, 0.59898, 0.922271, 1.12165, 1.22933,
        1.2831, 1.30859, 1.32024, 1.32542, 1.32768]) < 1e-5
end

function test_na()
  s = pgbf(1.0)
  c = cgbf(0.0,0.0,0.0)
  push!(c,1,1)
  @assert isapprox(amplitude(c,0,0,0),0.71270547)
  @assert isapprox(nuclear_attraction(s,s,0.,0.,0.),-1.59576912)
  @assert isapprox(nuclear_attraction(c,c,0.,0.,0.),-1.595769)
end

function test_fgamma()
  for (x,res) in Any[(0.,1),
          (30.,0.161802),
          (60.,0.1144114),
          (90.,0.0934165),
          (120.,0.08090108),
          (300.,0.051166336)]
    #@assert isapprox(res,Fgamma(0,x))
    @printf("res, Fgamma(0,x) %18.10f %18.10f\n", x, Fgamma(0,x))
  end
  @printf "test_fgamma is passed\n"
end

#todo make into a test
function test_one()
  s1 = pgbf(1)
  s2 = pgbf(1,0,1,0)
  x=y=0.
  println("S: $(overlap(s1,s2))")
  println("T: $(kinetic(s1,s2))")
  for z in linspace(0,1,5)
    println("V: $z $(nuclear_attraction(s1,s2,x,y,z))")
  end
end

function test_na2()
  li,h = lih.atomlist
  bfs = build_basis(lih)
  s1,s2,x,y,z,h1s = bfs.bfs
  @assert isapprox(nuclear_attraction(s1,s1,lih),-8.307532656)
end

nuclear_repulsion(a::Atom,b::Atom)= a.atno*b.atno/sqrt(dist2(a.x-b.x,a.y-b.y,a.z-b.z))
function nuclear_repulsion(mol::Molecule)
  nr = 0
  for (i,j) in pairs(nat(mol),"subdiag")
    nr += nuclear_repulsion(mol.atomlist[i],mol.atomlist[j])
  end
  return nr
end

function test_geo_basis()

  @printf("\nCalling test_geo_basis\n")

  E1     = nuclear_repulsion(h2)
  E1_ref = 0.7223600367
  @printf("H2: nuclear_repulsion %18.10f, err = %18.10f\n", E1, abs(E1_ref-E1))

  @assert nel(h2) == 2
  @assert nel(h2o) == 10
  @assert length(sto3g)==10

  bfs = build_basis(h2)
  @assert length(bfs.bfs)==2

  l,r = bfs.bfs
  @assert isapprox(overlap(l,l),1)
  @assert isapprox(overlap(r,r),1)

  #@assert isapprox(overlap(l,r),0.66473625)
  ovl1 = overlap(l,r)
  ovl1_ref = abs(ovl1 - 0.66473625)
  @printf("overlap(l,r), err %18.10f, %18.10f\n", ovl1, ovl1_ref)

  @assert isapprox(kinetic(l,l),0.76003188)
  @assert isapprox(kinetic(r,r),0.76003188)
  @assert isapprox(kinetic(l,r),0.24141861181119084)
  @assert isapprox(coulomb(l,l,l,l), 0.7746059439196398)
  @assert isapprox(coulomb(r,r,r,r), 0.7746059439196398)
  @assert isapprox(coulomb(l,l,r,r), 0.5727937653511646)
  @assert isapprox(coulomb(l,l,l,r), 0.4488373301593464)
  @assert isapprox(coulomb(l,r,l,r), 0.3025451156654606)
  bfs = build_basis(h2o)

  s1,s2,px,py,pz,hl,hr = bfs.bfs
  @assert isapprox(coulomb(s1,s2,hl,hr),0.03855344493645537)
  @assert isapprox(coulomb(s1,pz,hl,hr),-0.0027720110485359053)
  @assert isapprox(coulomb(s1,hl,pz,hr),-0.010049491284827426)
  @assert coulomb(s1,py,hl,hr)==0
  @assert coulomb(s1,hl,py,hr)==0
end


text:
# layers.jl

struct GCNLayer
    W
    b
end

GCNLayer(in::Integer, out::Integer) =
  GCNLayer(param(randn(out, in)), param(randn(out)))

(m::GCNLayer)(x, adj) = adj * m.W * x .+ m.b

max_stars_repo_name: mcbaron/fluxGCN

max_issues_repo_name: mcbaron/fluxGCN

max_forks_repo_name: mcbaron/fluxGCN


text:
module JungleHelperPredatorPlant

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/PredatorPlant" PredatorPlant(x::Integer, y::Integer, color::String="Pink", facingRight::Bool=false, sprite::String="")

const colors = String["Pink", "Blue", "Yellow"]

const placements = Ahorn.PlacementDict(
    "Predator Plant ($(color)) (Jungle Helper)" => Ahorn.EntityPlacement(
        PredatorPlant,
        "point",
        Dict{String, Any}(
          "color" => color
        )
    ) for color in colors
)

const spritePaths = Dict{String,String}(
    "Pink" => "JungleHelper/PredatorPlant/Bitey00",
    "Blue" => "JungleHelper/PredatorPlant/Blue/BiteyB00",
    "Yellow" => "JungleHelper/PredatorPlant/Yellow/BiteyC00"
)

Ahorn.editingOptions(entity::PredatorPlant) = Dict{String, Any}(
    "color" => colors
)

function Ahorn.selection(entity::PredatorPlant)
    x, y = Ahorn.position(entity)
    facingRight = get(entity.data, "facingRight", false)
    return Ahorn.Rectangle(x - (facingRight ? 18 : 6), y - 16, 24, 24)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::PredatorPlant, room::Maple.Room)
    color = get(entity.data, "color", "Pink")
    scaleX = get(entity.data, "facingRight", false) ? -1 : 1
    Ahorn.drawSprite(ctx, spritePaths[color], 0, -4, sx=scaleX)
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
using SimpleSolvers
using SimpleSolvers: LinesearchState, StaticState
using Test

include("optimizers_problems.jl")

f(x) = x^2-1
g(x) = 2x

function F(x)
    # sum(x.^2) - 1
    y = - one(eltype(x))
    for _x in x
        y += _x^2
    end
    return y
end

@testset "$(rpad("Bracketing",80))" begin
    @test bracket_minimum(x -> x^2) == (-SimpleSolvers.DEFAULT_BRACKETING_s, +SimpleSolvers.DEFAULT_BRACKETING_s)
    @test bracket_minimum(x -> (x-1)^2) == (0.64, 2.56)
end

@testset "$(rpad("Static",80))" begin
    # xᵤ = 3.
    x₀ = [-3.]
    x₁ = [+1.]
    δx = x₁ - x₀
    x  = copy(x₀)

    o  = MultivariateObjective(f, x)
    ls = StaticState()

    @test ls == LinesearchState(Static(), o)
    @test ls == LinesearchState(Static(), f, x)
    @test ls == LinesearchState(Static(), f, x; D = g)
    @test ls == LinesearchState(Static(1.0), o)
    @test ls == LinesearchState(Static(1.0), f, x)
    @test ls == LinesearchState(Static(1.0), f, x; D = g)

    @test ls == StaticState(o)
    @test ls == StaticState(f, x)
    @test ls == StaticState(f, x; D = g)

    x1 = copy(x₀); x2 = copy(x₀); @test solve!(x1, δx, ls) == ls(x2, δx) == x₁

    o1  = MultivariateObjective(f, x)
    o2  = MultivariateObjective(f, g, x)
    ls1 = Linesearch(x, o1; algorithm = Static())
    ls2 = Linesearch(x, o2; algorithm = Static())
    ls3 = Linesearch(x, o1; algorithm = Static(1.0))

    # @test ls1 == Linesearch(x, F; algorithm = Static())
    # @test ls2 == Linesearch(x, F; algorithm = Static(), D = D)

    x1 = copy(x₀); x2 = copy(x₀); @test solve!(x1, δx, ls1) == ls1(x2, δx) == x₁
    x1 = copy(x₀); x2 = copy(x₀); @test solve!(x1, δx, ls2) == ls2(x2, δx) == x₁
    x1 = copy(x₀); x2 = copy(x₀); @test solve!(x1, δx, ls3) == ls3(x2, δx) == x₁

end

@testset "$(rpad("Bisection",80))" begin

    x₀ = -3.0
    x₁ = +0.5
    δx = x₁ .- x₀
    x  = copy(x₀)

    x1 = bisection(f, x₀, x₁; config = Options(x_abstol = zero(x)))
    x2 = bisection(f, x; config = Options(x_abstol = zero(x)))

    @test x1 ≈ -1  atol=∛(2eps())
    @test x2 ≈ -1  atol=∛(2eps())

    x1 = bisection(f, x₀, x₁; config = Options(f_abstol = zero(f(x))))
    x2 = bisection(f, x; config = Options(f_abstol = zero(f(x))))

    @test x1 ≈ -1  atol=2eps()
    @test x2 ≈ -1  atol=2eps()

    # x₀ = [-3.0]
    # x₁ = [+0.5]
    # xₛ = [-1.0]
    # δx = x₁ .- x₀
    # x  = copy(x₀)

    # x_abstol = zero(eltype(x))
    # f_abstol = zero(F(x))

    # ls = Linesearch(x, F; algorithm = Bisection(), config = Options(x_abstol = x_abstol))

    # x1 = copy(x₀)
    # x2 = copy(x₀)
    # x3 = copy(x₀)

    # ls(x1, δx)
    # solve!(x2, δx, ls)
    # solve!(x3, δx, ls.state)

    # @test x1 ≈ xₛ  atol=∛(2eps())
    # @test x2 ≈ xₛ  atol=∛(2eps())
    # @test x3 ≈ xₛ  atol=∛(2eps())

    # ls = Linesearch(x, F; algorithm = Bisection(), config = Options(f_abstol = f_abstol))

    # x1 = copy(x₀)
    # x2 = copy(x₀)
    # x3 = copy(x₀)

    # ls(x1, δx)
    # solve!(x2, δx, ls)
    # solve!(x3, δx, ls.state)

    # @test x1 ≈ xₛ  atol=2eps()
    # @test x2 ≈ xₛ  atol=2eps()
    # @test x3 ≈ xₛ  atol=2eps()

end

@testset "$(rpad("Backtracking",80))" begin

    x₀ = -3.0
    x₁ = +3.0
    xₛ =  0.0
    δx = x₁ .- x₀

    _f = α -> f(x₀ + α * δx)
    _d = α -> g(x₀ + α * δx)

    ls = Linesearch(x₀, _f; D = _d, algorithm = Backtracking(), config = Options(x_abstol = zero(x₀)))

    α1 = ls()
    # solve!(x2, δx, ls)
    α2 = backtracking(_f, x₀; config = Options(x_abstol = zero(x₀)))
    α3 = backtracking(_f, x₀; D = _d, config = Options(x_abstol = zero(x₀)))

    x1 = x₀ + α1 * δx
    x2 = x₀ + α2 * δx
    x3 = x₀ + α3 * δx

    @test x1 ≈ xₛ  atol=∛(2eps())
    @test x2 ≈ xₛ  atol=∛(2eps())
    @test x3 ≈ xₛ  atol=∛(2eps())

    # x₀ = [-3.0]
    # x₁ = [+3.0]
    # xₛ = [ 0.0]
    # δx = x₁ .- x₀
    # x  = copy(x₀)

    # x1 = copy(x₀)
    # x2 = copy(x₀)
    # x3 = copy(x₀)

    # ls = Linesearch(x, F; algorithm = Backtracking(), config = Options(x_abstol = zero(eltype(x))))

    # ls(x1, δx)
    # solve!(x2, δx, ls)
    # backtracking(F, x3, δx; config = Options(x_abstol = zero(eltype(x))))

    # @test x1 ≈ xₛ  atol=∛(2eps())
    # @test x2 ≈ xₛ  atol=∛(2eps())
    # @test x3 ≈ xₛ  atol=∛(2eps())

    # n = 1
    # x₀ = -0.5*ones(n)
    # x₁ = +1.0*ones(n)
    # x  = copy(x₀)
    # f  = zeros(n)
    # g  = zeros(n,n)
    # ls = Solver(F!, x, f)

    # F!(f,x)
    # J!(g,x)

    # solve!(x, f, g, x₀, x₁, ls)

    # F!(f,x)

    # @test x ≈ zero(x) atol=4E-1
    # @test f ≈ zero(f) atol=8E-2

    # @test solve!(x, f, g, x₀, x₁, ls) == ls(x, f, g, x₀, x₁) == SFunc(F!, x, f, g, x₀, x₁)

    # n = 3
    # x₀ = -ones(n)
    # x₁ = +ones(n)
    # x  = copy(x₀)
    # f  = zeros(n)
    # g  = zeros(n,n)
    # ls = Solver(F!, x, f)

    # F!(f,x)
    # J!(g,x)

    # solve!(x, f, g, x₀, x₁, ls)

    # F!(f,x)

    # @test x == zero(x)
    # @test f == zero(f)

    # @test solve!(x, f, g, x₀, x₁, ls) == ls(x, f, g, x₀, x₁) == SFunc(F!, x, f, g, x₀, x₁)

end

@testset "$(rpad("Polynomial Linesearch",80))" begin

    x₀ = -3.0
    x₁ = +3.0
    xₛ =  0.0
    δx = x₁ .- x₀

    _f = α -> f(x₀ + α * δx)
    _d = α -> g(x₀ + α * δx)

    ls = Linesearch(x₀, _f; D = _d, algorithm = Quadratic(), config = Options(x_abstol = zero(x₀)), σ₀=0.5, σ₁=0.8)

    α1 = ls()
    # solve!(x2, δx, ls)
    α2 = quadratic(_f, x₀; config = Options(x_abstol = zero(x₀)), σ₀=0.5, σ₁=0.8)
    α3 = quadratic(_f, x₀; D = _d, config = Options(x_abstol = zero(x₀)), σ₀=0.5, σ₁=0.8)

    x1 = x₀ + α1 * δx
    x2 = x₀ + α2 * δx
    x3 = x₀ + α3 * δx

    @test x1 ≈ xₛ  atol=∛(2eps())
    @test x2 ≈ xₛ  atol=∛(2eps())
    @test x3 ≈ xₛ  atol=∛(2eps())

    # x₀ = [-3.0]
    # x₁ = [+3.0]
    # xₛ = [ 0.0]
    # δx = x₁ .- x₀
    # x  = copy(x₀)

    # x1 = copy(x₀)
    # x2 = copy(x₀)
    # x3 = copy(x₀)

    # ls = Linesearch(x, F; algorithm = Quadratic(), config = Options(x_abstol = zero(eltype(x))))

    # ls(x1, δx)
    # solve!(x2, δx, ls)
    # quadratic(F, x3, δx; config = Options(x_abstol = zero(eltype(x))))

    # @test x1 ≈ xₛ  atol=∛(2eps())
    # @test x2 ≈ xₛ  atol=∛(2eps())
    # @test x3 ≈ xₛ  atol=∛(2eps())

    # n = 1
    # x₀ = -0.5*ones(n)
    # x₁ = +1.0*ones(n)
    # x  = copy(x₀)
    # f  = zeros(n)
    # g  = zeros(n,n)
    # ls = Solver(F!, x, f)

    # F!(f,x)
    # J!(g,x)

    # solve!(x, f, g, x₀, x₁, ls)

    # F!(f,x)

    # @test x ≈ zero(x) atol=4E-1
    # @test f ≈ zero(f) atol=8E-2

    # @test solve!(x, f, g, x₀, x₁, ls) == ls(x, f, g, x₀, x₁) == SFunc(F!, x, f, g, x₀, x₁)

    # n = 3
    # x₀ = -ones(n)
    # x₁ = +ones(n)
    # x  = copy(x₀)
    # f  = zeros(n)
    # g  = zeros(n,n)
    # ls = Solver(F!, x, f)

    # F!(f,x)
    # J!(g,x)

    # solve!(x, f, g, x₀, x₁, ls)

    # F!(f,x)

    # @test x == zero(x)
    # @test f == zero(f)

    # @test solve!(x, f, g, x₀, x₁, ls) == ls(x, f, g, x₀, x₁) == SFunc(F!, x, f, g, x₀, x₁)

end


text:
"""
    internal_rate_of_return(cashflows::vector)
    internal_rate_of_return(cashflows::Vector, timepoints::Vector)

Calculate the internal_rate_of_return with given timepoints. If no timepoints given, will assume that a series of equally spaced cashflows, assuming the first cashflow occurring at time zero.

Will try to return a root within the range [-2,2]. If the fast solver does not find one matching this condition, then a more robust search will be performed over the [.99,2] range.

The solution returned will be in the range [-2,2], but may not be the one nearest zero. For a slightly slower, but more robust version, call `ActuaryUtilities.irr_robust(cashflows,timepoints)` directly.

# Example
```julia-repl
julia> internal_rate_of_return([-100,110],[0,1]) # e.g. cashflows at time 0 and 1
0.10000000001652906
julia> internal_rate_of_return([-100,110]) # implied the same as above
0.10000000001652906
```

"""
function internal_rate_of_return(cashflows)

    return internal_rate_of_return(cashflows, 0:length(cashflows)-1)

end

function internal_rate_of_return(cashflows,times)
    # first try to quickly solve with newton's method, otherwise
    # revert to a more robust method
    lower,upper = -2.,2.

    v = try
        irr_newton(cashflows,times)
    catch e
        if isa(e,Roots.ConvergenceFailed) || sprint(showerror, e) =="No convergence"
            return irr_robust(cashflows,times)
        else
            throw(e)
        end
    end

    if v <= upper && v >= lower
        return v
    else
        return irr_robust(cashflows,times)
    end
end

irr_robust(cashflows) = irr_robust(cashflows,0:length(cashflows)-1)

function irr_robust(cashflows, times)
    f(i) =  sum(cf / (1+i)^t for (cf,t) in zip(cashflows,times))
    # lower bound at -.99 because otherwise we can start taking the root of a negative number
    # when a time is fractional.
    roots = Roots.find_zeros(f, -0.99, 2)

    # short circuit and return nothing if no roots found
    isempty(roots) && return nothing
    # find and return the one nearest zero
    min_i = argmin(roots)
    return roots[min_i]

end

irr_newton(cashflows) = irr_newton(cashflows,0:length(cashflows)-1)

function irr_newton(cashflows, times)
    # use newton's method with hand-coded derivative
    f(r) =  sum(cf * exp(-r*t) for (cf,t) in zip(cashflows,times))
    f′(r) = sum(-t*cf * exp(-r*t) for (cf,t) in zip(cashflows,times) if t > 0)
    # r = Roots.solve(Roots.ZeroProblem((f,f′), 0.0), Roots.Newton())
    r = Roots.newton(x->(f(x),f(x)/f′(x)),0.0)
    return exp(r)-1

end

"""
    irr()

    An alias for `internal_rate_of_return`.
"""
irr = internal_rate_of_return

"""
    present_value(interest, cashflows::Vector, timepoints)
    present_value(interest, cashflows::Vector)

Discount the `cashflows` vector at the given `interest_interestrate`,  with the cashflows occurring
at the times specified in `timepoints`. If no `timepoints` given, assumes that cashflows happen at times 1,2,...,n.

The `interest` can be an `InterestCurve`, a single scalar, or a vector wrapped in an `InterestCurve`.

# Examples
```julia-repl
julia> present_value(0.1, [10,20],[0,1])
28.18181818181818
julia> present_value(Yields.Forward([0.1,0.2]), [10,20],[0,1])
28.18181818181818 # same as above, because first cashflow is at time zero
```

Example on how to use real dates using the [DayCounts.jl](https://github.com/JuliaFinance/DayCounts.jl) package
```jldoctest

using DayCounts
dates = Date(2012,12,31):Year(1):Date(2013,12,31)
times = map(d -> yearfrac(dates[1], d, DayCounts.Actual365Fixed()),dates) # [0.0,1.0]
present_value(0.1, [10,20],times)

# output
28.18181818181818

```

"""
function present_value(yc::T, cashflows, timepoints) where {T <: Yields.AbstractYield}
    sum(discount(yc,t) * cf for (t,cf) in zip(timepoints, cashflows))
end

function present_value(yc::T, cashflows) where {T <: Yields.AbstractYield}
    present_value(yc,cashflows,1:length(cashflows))
end

function present_value(i, x)

    v = 1.0
    v_factor = 1/(1+i)
    pv = 0.0

    for (t,cf) in zip(1:length(x),x)
        v *= v_factor
        pv += v * cf
    end
    return pv
end

function present_value(i, v, times)
    return present_value(Yields.Constant(i), v, times)
end

# Interest Given is an array, assume forwards.
function present_value(i::AbstractArray, v)
    yc = Yields.Forward(i)
    return sum(discount(yc, t) * cf for (t,cf) in zip(1:length(v),v))
end

# Interest Given is an array, assume forwards.
function present_value(i::AbstractArray, v, times)
    yc = Yields.Forward(i, times)
    return sum(discount(yc, t) * cf for (cf, t) in zip(v,times))
end

"""
    pv()

    An alias for `present_value`.
"""
pv = present_value

"""
    present_value(interest, cashflows::Vector, timepoints)
    present_value(interest, cashflows::Vector)

Efficiently calculate a vector representing the present value of the given cashflows at each period prior to the given timepoint.

# Examples
```julia-repl
julia> present_values(0.00, [1,1,1])
[3,2,1]

julia> present_values(Yields.Forward([0.1,0.2]), [10,20],[0,1])
2-element Vector{Float64}:
 28.18181818181818
 18.18181818181818
```

"""
function present_values(interest, cashflows)
    pvs = Vector{Float64}(undef,length(cashflows))
    pvs[end] = Yields.discount(interest, lastindex(cashflows) - 1, lastindex(cashflows)) * cashflows[end]
    for (t, cf) in Iterators.reverse(enumerate(cashflows[1:end - 1]))
        pvs[t] = Yields.discount(interest, t - 1, t) * (cf + pvs[t + 1])
    end

    return pvs
end

function present_values(interest,cashflows,times)
    present_values_accumulator(interest,cashflows,times)
end

function present_values_accumulator(interest,cashflows,times,pvs=[0.0])
    from_time = length(times) == 1 ? 0. : times[end-1]
    pv = discount(interest,from_time,last(times)) *(first(pvs) + last(cashflows))
    pvs = pushfirst!(pvs,pv)

    if length(cashflows) > 1

        new_cfs = @view cashflows[1:end-1]
        new_times = @view times[1:end-1]
        return present_values_accumulator(interest,new_cfs,new_times,pvs)
    else
        # last discount and return
        return pvs[1:end-1] # end-1 get rid of trailing 0.0
    end
end

# if given a vector of rates, assume that it should be a forward discount yield
function present_values(y::Vector{T}, cfs, times) where {T <: Real}
    return present_values(Yields.Forward(y), cfs, times)
end

"""
    price(...)

The absolute value of the `present_value(...)`.

# Extended help

Using `price` can be helpful if the directionality of the value doesn't matter. For example, in the common usage, duration is more interested in the change in price than present value, so `price` is used there.
"""
price(x1,x2) = present_value(x1, x2) |> abs
price(x1,x2,x3) = present_value(x1, x2, x3) |> abs

"""
    breakeven(yield, cashflows::Vector)
    breakeven(yield, cashflows::Vector,times::Vector)

Calculate the time when the accumulated cashflows breakeven given the yield.

Assumptions:

- cashflows occur at the end of the period
- cashflows evenly spaced with the first one occuring at time zero if `times` not given

Returns `nothing` if cashflow stream never breaks even.

```jldoctest
julia> breakeven(0.10, [-10,1,2,3,4,8])
5

julia> breakeven(0.10, [-10,15,2,3,4,8])
1

julia> breakeven(0.10, [-10,-15,2,3,4,8]) # returns the `nothing` value

```
"""
function breakeven(y::T, cashflows::Vector, timepoints::Vector) where {T <: Yields.AbstractYield}
    accum = zero(eltype(cashflows))
    last_neg = nothing

    accum += cashflows[1]
    if accum >= 0 && isnothing(last_neg)
        last_neg = timepoints[1]
    end

    for i in 2:length(cashflows)
        # accumulate the flow from each timepoint to the next
        accum *= Yields.accumulation(y, timepoints[i - 1], timepoints[i])
        accum += cashflows[i]

        if accum >= 0 && isnothing(last_neg)
            last_neg = timepoints[i]
        elseif accum < 0
            last_neg = nothing
        end
    end

    return last_neg

end

function breakeven(y::T, cfs, times) where {T <: Real}
    return breakeven(Yields.Constant(y), cfs, times)
end

function breakeven(y::Vector{T}, cfs, times) where {T <: Real}
    return breakeven(Yields.Forward(y), cfs, times)
end

function breakeven(i, cashflows::Vector)
    return breakeven(i, cashflows, [t for t in 0:length(cashflows) - 1])
end

abstract type Duration end

struct Macaulay <: Duration end
struct Modified <: Duration end
struct DV01 <: Duration end

"""
    duration(Macaulay(),interest_rate,cfs,times)
    duration(Modified(),interest_rate,cfs,times)
    duration(::DV01,interest_rate,cfs,times)
    duration(interest_rate,cfs,times)             # Modified Duration
    duration(interest_rate,valuation_function)    # modified Duration

Calculates the Macaulay, Modified, or DV01 duration. `times` may be ommitted and the valuation will assume evenly spaced cashflows starting at the end of the first period.
- `interest_rate` should be a fixed effective yield (e.g. `0.05`).

When not given `Modified()` or `Macaulay()` as an argument, will default to `Modified()`.

# Examples

Using vectors of cashflows and times
```julia-repl
julia> times = 1:5
julia> cfs = [0,0,0,0,100]
julia> duration(0.03,cfs,times)
4.854368932038834
julia> duration(Macaulay(),0.03,cfs,times)
5.0
julia> duration(Modified(),0.03,cfs,times)
4.854368932038835
julia> convexity(0.03,cfs,times)
28.277877274012614

```

Using any given value function:

```julia-repl
julia> lump_sum_value(amount,years,i) = amount / (1 + i ) ^ years
julia> my_lump_sum_value(i) = lump_sum_value(100,5,i)
julia> duration(0.03,my_lump_sum_value)
4.854368932038835
julia> convexity(0.03,my_lump_sum_value)
28.277877274012617

```
"""
function duration(::Macaulay, yield, cfs, times)
    return sum(times .* price.(yield, vec(cfs), times) / price(yield, vec(cfs), times))
end

function duration(::Modified, yield, cfs, times)
    D(i) = price(i, vec(cfs), times)
    return duration(yield, D)
end

function duration(yield, valuation_function)
    D(i) = log(valuation_function(i + yield))
    δV =  - ForwardDiff.derivative(D, 0.0)
end

function duration(yield::Y, valuation_function) where {Y <: Yields.AbstractYield}
    D(i) = log(valuation_function(i + yield))
    δV =  - ForwardDiff.derivative(D, 0.0)
end

function duration(yield, cfs, times)
    return duration(Modified(), yield, vec(cfs), times)
end
function duration(yield::Y, cfs::A) where {Y <: Yields.AbstractYield,A <: AbstractArray}
    times = 1:length(cfs)
    return duration(Modified(), yield, vec(cfs), times)
end

function duration(yield::R, cfs) where {R <: Real}
    return duration(Yields.Constant(yield), cfs)
end

function duration(::DV01, yield, cfs, times)
    return duration(DV01(), yield, i -> price(i, vec(cfs), times))
end
function duration(d::Duration, yield, cfs)
    times = 1:length(cfs)
    return duration(d, yield, vec(cfs), times)
end

function duration(::DV01, yield, valuation_function)
    return duration(yield, valuation_function) * valuation_function(yield) / 100
end

"""
    convexity(yield,cfs,times)
    convexity(yield,valuation_function)

Calculates the convexity.
    - `yield` should be a fixed effective yield (e.g. `0.05`).
    - `times` may be omitted and it will assume `cfs` are evenly spaced beginning at the end of the first period.

# Examples

Using vectors of cashflows and times
```julia-repl
julia> times = 1:5
julia> cfs = [0,0,0,0,100]
julia> duration(0.03,cfs,times)
4.854368932038834
julia> duration(Macaulay(),0.03,cfs,times)
5.0
julia> duration(Modified(),0.03,cfs,times)
4.854368932038835
julia> convexity(0.03,cfs,times)
28.277877274012614

```

Using any given value function:

```julia-repl
julia> lump_sum_value(amount,years,i) = amount / (1 + i ) ^ years
julia> my_lump_sum_value(i) = lump_sum_value(100,5,i)
julia> duration(0.03,my_lump_sum_value)
4.854368932038835
julia> convexity(0.03,my_lump_sum_value)
28.277877274012617

```

"""
function convexity(yield, cfs, times)
    return convexity(yield, i -> price(i, vec(cfs), times))
end

function convexity(yield, cfs::A) where {A <: AbstractArray}
    times = 1:length(cfs)
    return convexity(yield, i -> price(i, vec(cfs), times))
end

function convexity(yield, valuation_function)
    v(x) = abs(valuation_function(yield + x[1]))
    ∂²P = ForwardDiff.hessian(v, [0.0])
    return ∂²P[1] / v([0.0])
end


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 1.0.5
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# ### 2D+T dynamic MRI simulation
# using GA radial sampling
# reconstructed with temporal "TV" regularizer (corner-rounded)
# 2019-06-13, Jeff Fessler
# 2019-06-23 update to use more realistic simulated sensitivity maps
# 2020-06-20 update

using MIRT: jim, prompt
using MIRT: image_geom, ellipse_im_params, ellipse_im
using MIRT: nufft_init, diffl_map, ncg
using MIRT: ir_mri_sensemap_sim, ir_mri_kspace_ga_radial
using Plots: gui, plot, scatter, default; default(markerstrokecolor=:auto)
using LinearAlgebra: norm, dot, Diagonal
using LinearMapsAA
using Random: seed!
jim(:abswarn, false); # suppress warnings about display of |complex| images

# generate dynamic image sequence
if !@isdefined(xtrue)
	N = (60,64)
	fov = 220
	nt = 8 # frames
	ig = image_geom(nx=N[1], ny=N[2], fov=fov)

	ellpar = ellipse_im_params(ig, :southpark)
	ellpars = Array{Float32}(undef, size(ellpar)..., nt)
	xtrue = Array{ComplexF32}(undef, N..., nt)
	for it=1:nt
		tmp = ellpar
		tmp[2,4] = 15 + 5 * sin(2*pi*it/nt) # mouth open/close
		ellpars[:,:,it] = tmp
		xtrue[:,:,it] = ellipse_im(ig, ellpars[:,:,it], oversample=4)
		jim(ig.x, ig.y, xtrue[:,:,it], title="frame $it")
		gui()
	end
end
jim(xtrue, yflip=ig.dy < 0)
#prompt()

# plot one time course to see temporal change
ix,iy = 30,14
plot(1:nt, abs.(xtrue[ix,iy,:]), label="ix=$ix, iy=$iy",
	marker=:o, xlabel="frame")
#prompt()

# +
# k-space sampling and data
if !@isdefined(kspace)
	@show accelerate = 3
	@show nspf = round(Int, maximum(N)/accelerate) # spokes per frame
	Nro = maximum(N)
	Nspoke = nspf * nt
	kspace = ir_mri_kspace_ga_radial(Nro = Nro, Nspoke = Nspoke)
	kspace[:,:,1] ./= ig.fovs[1]
	kspace[:,:,2] ./= ig.fovs[2]
	kspace = reshape(kspace, Nro, nspf, nt, 2)

	if true # plot sampling (in units of cycles/pixel)
		ps = Array{Any}(undef, nt)
		for it=1:nt
			ps[it] = scatter(kspace[:,:,it,1] * ig.fovs[1], kspace[:,:,it,2] * ig.fovs[2],
				xtick=(-1:1)*0.5, ytick=(-1:1)*0.5, xlim=[-1,1]*0.52, ylim=[-1,1]*0.52,
				aspect_ratio=1, label="", markersize=1)
			plot(ps[it])
			gui()
		end
		plot(ps..., layout=(2,4))
	#	prompt()
	end
end
# -

#+
# make sensitivity maps, normalized so SSoS = 1
if !@isdefined(smap)
	ncoil = 2
	smap = ir_mri_sensemap_sim(dims=N, ncoil=ncoil, orbit_start=[90])
	p1 = jim(smap, "ncoil=$ncoil sensitivity maps raw")

	ssos = sqrt.(sum(abs.(smap).^2, dims=ndims(smap))) # SSoS
	ssos = selectdim(ssos, ndims(smap), 1)
	p2 = jim(ssos, "SSoS")

	for ic=1:ncoil
		selectdim(smap, ndims(smap), ic) ./= ssos
	end
	p3 = jim(smap, "ncoil=$ncoil sensitivity maps")

	ssos = sqrt.(sum(abs.(smap).^2, dims=ndims(smap))) # SSoS
	@assert all(isapprox.(ssos,1))
	plot(p1, p2, p3)
#	prompt()
end
#-

#+
# make system matrix for dynamic non-Cartesian parallel MRI
if !@isdefined(A)
	# a NUFFT object for each frame
	ns = Array{Any}(undef, nt)
	for it=1:nt
		om = [kspace[:,:,it,1][:] kspace[:,:,it,2][:]] * fov * 2 * pi
		ns[it] = nufft_init(om, N, n_shift = collect(N)/2)
	end

	# block diagonal system matrix, with one NUFFT per frame
	S = [Diagonal(selectdim(smap, ndims(smap), ic)[:]) for ic=1:ncoil]
	SO = s -> LinearMapAA(s ; idim=N, odim=N) # LinearMapAO
	AS1 = A1 -> vcat([A1 * SO(s) for s in S]...) # [A1*S1; ... ; A1*Sncoil]

	# output is essentially [nt Ncoil nspf Nro] (which is possibly unusual)
	# input is [N... nt]
	A = block_diag([AS1(s.A) for s in ns]...)
end
#-

#+
# simulate k-space data via an inverse crime
if !@isdefined(y)
	ytrue = A * xtrue

	snr2sigma = (db, yb) -> # compute noise sigma from SNR (no sqrt(2) needed)
		10^(-db/20) * norm(yb) / sqrt(length(yb))

	sig = Float32(snr2sigma(50, ytrue))
	seed!(0)
	y = ytrue + sig * randn(ComplexF32, size(ytrue))
	@show 20*log10(norm(ytrue) / norm(y - ytrue)) # verify SNR
end
#-

# initial image via zero-fill and scaling
if !@isdefined(x0)
	# todo: should use density compensation, perhaps via
	# https://github.com/JuliaGeometry/VoronoiDelaunay.jl
	x0 = A' * y # zero-filled recon (for each frame)
	tmp = A * x0 # Nkspace × Ncoil × Nframe
	x0 = (dot(tmp,y) / norm(tmp)^2) * x0 # scale sensibly
	jim(x0, "initial image", yflip=ig.dy < 0)
end
#prompt()

# +
# temporal finite differences
if !@isdefined(Dt)
	Dt = diffl_map((N..., nt), length(N)+1 ; T=eltype(A))
	tmp = Dt' * (Dt * xtrue)
	jim(tmp, "time diff", yflip=ig.dy < 0)
end
#prompt()
# -

# +
# run nonlinear CG
if !@isdefined(xh)
	niter = 90
	delta = Float32(0.1) # small relative to temporal differences
	reg = Float32(2^20) # trial and error here
	ffair = (t,d) -> d^2 * (abs(t)/d - log(1 + abs(t)/d))
	pot = z -> ffair(z, delta)
	dpot = z -> z / (Float32(1) + abs(z/delta))
	cost = x -> 0.5 * norm(A*x - y)^2 + reg * sum(pot.(Dt * x))
	fun = (x,iter) -> cost(x)
	gradf = [v -> v - y, u -> reg * dpot.(u)]
	curvf = [v -> Float32(1), u -> reg]
	B = [A, Dt]
	(xh, out) = ncg(B, gradf, curvf, x0 ; niter=niter, fun=fun)
	costs = [out[i+1][1] for i=0:niter]
end

# show results
	plot(layout=(2,2),
		jim(xtrue, "xtrue", yflip=ig.dy < 0),
		jim(xh, "recon", yflip=ig.dy < 0),
		jim(xh-xtrue, "error", yflip=ig.dy < 0),
		scatter(0:niter, log.(costs), label="cost", xlabel="iteration"),
	)

max_stars_repo_name: cecroc/mirt-demo

max_issues_repo_name: cecroc/mirt-demo

max_forks_repo_name: cecroc/mirt-demo


text:
using UnregisteredWithoutProject
using Test

@test f(2) == 2

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
# Make sure all solver options can be created without specifying float type
@test_nowarn ilqr = iLQRSolverOptions()
@test_nowarn al = AugmentedLagrangianSolverOptions()
@test_nowarn altro = ALTROSolverOptions()
@test_nowarn pn = ProjectedNewtonSolverOptions()


text:
#addprocs(Sys.CPU_CORES)
#@everywhere using Plasm

using Plasm
using PyCall
p = PyCall.pyimport("pyplasm")

#####

ncubes = 3
V,cells = LinearAlgebraicRepresentation.larCuboids([ncubes,ncubes,ncubes],true)
VV,EV,FV,CV = cells
Plasm.viewexploded(V,FV);

t = -ncubes/2
V = Plasm.translate([t,t,t],V)
Plasm.viewexploded(V,FV);

W = copy(V)
FW = copy(FV)

W = Plasm.rotate((0,π/3,0),Plasm.rotate((π/3,0,0), W))
Plasm.viewexploded(W,FW);

V,W = V',W'
EW = characteristicMatrix(EV)
FE = boundary2(FV,EV)

rubik = [V,EW,FE]
rot_rubik = [W,EW,FE]
two_rubiks = LinearAlgebraicRepresentation.skel_merge(rubik..., rot_rubik...)
arranged_rubiks = LinearAlgebraicRepresentation.spatial_arrangement(two_rubiks...,multiproc=false)

V,cscEV,cscFE,cscCF = arranged_rubiks

ne,nv = size(cscEV)
EV = [findn(cscEV[e,:]) for e=1:ne]
Plasm.viewexploded(V',EV)

nf = size(cscFE,1)
FV = [collect(Set(vcat([EV[e] for e in findn(cscFE[f,:])]...)))  for f=1:nf]
Plasm.viewexploded(V',FV)

nc = size(cscCF,1)
CV = [collect(Set(vcat([FV[f] for f in findn(cscCF[c,:])]...)))  for c=2:nc]
Plasm.viewexploded(V',CV)

######

V,(VV,EV,FV,CV) = LinearAlgebraicRepresentation.larCuboids([2,2,1],true)
V,bases,coboundaries = LinearAlgebraicRepresentation.chaincomplex(V,FV,EV)
EV,FV,CV = bases
cscEV,cscFE,cscCF = coboundaries
Plasm.viewexploded(V,EV)
Plasm.viewexploded(V,FV)
Plasm.viewexploded(V,CV)
(ne,nv),nf,nc = size(cscEV),size(cscFE,1),size(cscCF,1)
nv-ne+nf-nc

#####

V,(VV,EV,FV,CV) = LinearAlgebraicRepresentation.larCuboids([2,2,1],true)
W,FW,EW = copy(V),copy(FV),copy(EV)
collection = [[W,FW,EW]]
for k=1:10
	W,FW,EW = rotate([0,0,π/15],copy(W)+.5),copy(FV),copy(EV)
	append!(collection, [[W,FV,EV]])
end
V,FV,EV = LinearAlgebraicRepresentation.collection2model(collection)
V,bases,coboundaries = LinearAlgebraicRepresentation.chaincomplex(V,FV,EV)
EV,FV,CV = bases
cscEV,cscFE,cscCF = coboundaries
Plasm.viewexploded(V,EV)
Plasm.viewexploded(V,FV)
Plasm.viewexploded(V,CV)

####

TV = LinearAlgebraicRepresentation.triangulate((1:length(FV),ones(length(FV))),V,FV,EV,cscFE,cscCF)
Plasm.viewexploded(V,TV)

########

Plasm.viewsolidcells(1.5,1.5,3.)(V,CV,FV,EV,cscCF,cscFE)


text:
# module multimode_steadystate

using QuantumOptics

function Halpha(system::MultimodeSystem, alpha::Number)
    H = Hamiltonian(system.particles)
    for n=1:length(system.modes)
        An = A(n, system)
        Bn = B(n, system)
        mode = system.modes[n]
        #H += 2*mode.eta*real(alpha)*Bn
        H += mode.eta^2*Bn^2*mode.delta/(mode.delta^2 + mode.kappa^2)
    end
    return H
end

function Hfunctional(system::MultimodeSystem, psi::Ket)
    E = expect(Hamiltonian(system.particles), psi)
    for n=1:length(system.modes)
        An = A(n, system)
        Bn = expect(B(n, system), psi)
        mode = system.modes[n]
        #H += 2*mode.eta*real(alpha)*Bn
        E += mode.eta^2*Bn^2*mode.delta/(mode.delta^2 + mode.kappa^2)
    end
    println("E = ", real(E))
    return real(E)
end

function Hξ(system::MultimodeSystem, ξ::Float64)
    H = Hamiltonian(system.particles)
    for n=1:length(system.modes)
        Bn = B(n, system)
        H += ξ*Bn
    end
    return H
end

# function steady_particle(system::MultimodeSystem, alpha<:Number...)
#     Hplus = Hamiltonian(system.particles)
#     Hminus = Hamiltonian(system.particles)
#     for n=1:length(system.modes)
#         An = A(n, system)
#         Bn = B(n, system)
#         mode = system.modes[n]
#         Hplus += mode.U0*abs2(alpha[n]*An) + 2*mode.eta*real(alpha[n])*Bn
#         Hminus += mode.U0*abs2(alpha[n]*An) - 2*mode.eta*real(alpha[n])*Bn
#     end

# end

# function steady_modes(system::MultimodeSystem, xplus::Ket)
#     alpha = Complex128[]
#     for n=1:length(system.modes)
#         a = dagger(xplus)*A(n, system)*xplus
#         b = dagger(xplus)*B(n, system)*xplus
#         push(alpha, b0)
#     end
# end

# end #multimode_steadystate


text:
#=
    A port of my Perl/Java/Python MakeRegex in Julia.

    make_regex(words)
    generates a regex for words. It use a simple approach which
    combined common prefixes in generating the regex.

    Some examples:
    * words = ["a", "al", "all", "alla", "an", "ann", "anna", "annas", "ananas"]
      regex: a(l(la?)?|n(anas|n(as?)?)?)?

    * words: ["and", "at", "do", "end", "for", "in", "is", "not", "of", "or", "use"]
      regex: (a(nd|t)|do|end|for|i[ns]|not|o[fr]|use)

    There is a simple way of handling character classes
    * words: ["price1", "price2", "price3", "price4"]
      regex: price[1234]

    If there is no common prefix then it just put '|' between the words
    * words: ["this", "is", "a", "very", "boring", "example", "with", "no", "common", "prefix"]
      regex: (a|boring|common|example|is|no|prefix|this|very|with)

    Also, see the (very old) page for my Perl package MakeRegex: http://hakank.org/makeregex/index.html
    The REAME file in that package states:
    """
    The Perl package MakeRegex composes a regex-expression from a list of
    words. It had been inspired by the emacs elisp module make-regex.el,
    by Simon Marshall.
    """

    This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
    See also my Julia page: http://www.hakank.org/julia/

=#

#
# common_prefix(p,  list)
#
# Here is where the main work is done. It's somewhat magically ported from my
# Perl/Java/Python versions...
#
function common_prefix(p, list)
    list_len = length(list)
    if list_len == 0
        return ""
    end
    if list_len == 1
        return p * join(list,"")
    end

    #
    # fix for some - as of now - unknown bug. To fix in later version!
    #
    if p == "" && list[1] == "" && list[2] == ""
        return ""
    end

    #
    # * Collect all the strings with the same prefix-char
    #
    hash = Dict()
    for word in sort(list)
        prefix = suffixed_word = ""
        if length(word) > 0
            prefix, suffixed_word... = word
        end
        # put the suffix in the list of other suffixes for
        # the this prefix
        hash[prefix] = push!(get(hash,prefix,[]),suffixed_word)
    end

    #
    # And recurse this list
    #
    all = []
    for key in keys(hash)
        comm = ""
        values = hash[key]
        if length(key) > 0
            sort!(values)
            comm = common_prefix(key, values)
        end

        # hack to be able to use the '?' char . Should be re-written!
        if comm == ""
            comm = " "
        end
        push!(all,comm)
    end
    sort!(all)

    # paren: what to put in parenthesis ('()' or '[]') if anything
    paren = ""
    all_len = length(all)
    if  all_len == 1
        paren = join(all,"")
    else
        len = maximum(length.(all))
        joinChar = len != 1 ? '|' : ""

        # joins all entries except for " "
        join_str = mark = ""
        count = 0
        for w in all
            got_hack_mark = w == " " ? true : false # This is a hack for handling '?'
            if length(w) > 0 && w != " "
                join_str *= w
                if count < all_len-1
                    join_str *=  joinChar
                end
            end

            if got_hack_mark
                mark = '?'
            end

            count = count + 1
        end

        paren = ""
        if length(join_str) === 1
            paren = join_str * mark
        else
            if len == 1
                paren = '[' * join_str * ']' * mark
            else
                paren = '(' * join_str * ')' * mark
            end
        end
    end
    return p * paren

end

function make_regex(words)
    replace.(words,r"([*?+])"=>s"\\\1") # replace meta characters
    # We sort the words to induce more common prefixes
    return common_prefix("", sort(words))

end

#
# check_regex(regex, words)
#
# Checks the regex againts a list of words.
#
function check_regex(regex, words)
    p = Regex(regex)
    for word in words
        println(word, " matches", !occursin(p,word) ? " NOT!" : "")
    end
end

tests = [
    ["all","alla"],

    # A lot of Swedish words
    [ "all", "alla", "alle", "alls", "palle", "palla", "pelle", "perkele",
      "ann", "anna", "annas", "anders", "håkan", "ångest", "ärlig",
      "solsken", "sture", "stina", "hörapparat", "hörsel", "hårig"],

    ["alla", "palla", "balla", "kalla", "all", "pall", "ball", "kall"],

    # "ananas" is the Swedish word for pineapple
    ["a", "al", "all", "alla", "an", "ann", "anna", "annas", "ananas"],

    ["a", "an", "ann", "anna", "annan", "annas", "annans", "ananas", "ananasens"],

    ["a", "ab", "abc", "abcd", "abcde", "abcdef", "b", "bc", "bcd", "bcde", "bcdef",
     "bcdefg", "abb", "abbc", "abbcc", "abbccdd"],

    ["this", "is", "a", "very", "boring", "example", "with", "no", "common", "prefix"],

    ["price1","price2","price3","price4"],

    # This is from Marshall's make-regex.el
    ["and", "at", "do", "end", "for", "in", "is", "not", "of", "or", "use"],

    # This is from Marshall's make-regex.el
    ["cond", "if", "while", "let*?", "prog1", "prog2", "progn",
     "catch", "throw", "save-restriction", "save-excursion",
     "save-window-excursion", "save-match-data", "unwind-protect",
     "condition-case", "track-mouse"],

    # This is from Marshall's make-regex.el
    ["abort", "abs", "accept", "access", "array",
      "begin", "body", "case", "constant", "declare",
      "delay", "delta", "digits", "else", "elsif", "entry",
      "exception", "exit", "function", "generic", "goto",
      "if", "others", "limited", "loop", "mod", "new",
      "null", "out", "subtype", "package", "pragma",
      "private", "procedure", "raise", "range", "record",
      "rem", "renames", "return", "reverse", "select",
      "separate", "task", "terminate", "then", "type",
      "when", "while", "with", "xor"]
]

for t in tests
    println("testing $t")
    println(make_regex(t))
    println()
end

#=
words = last(tests)
rx = make_regex(words)
println(words)
println("regex:$rx")
check_regex(rx,words)
=#

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# for sum(xs, weights) (#522)
Base.sum(xs::AbstractArray, weights::AbstractArray) = dot(xs, weights)
struct SumRuleConfig <: RuleConfig{Union{HasReverseMode}} end

const CFG = ChainRulesTestUtils.ADviaRuleConfig()

@testset "Reductions" begin
    @testset "sum(::Tuple)" begin
        test_frule(sum, Tuple(rand(5)))
    end
    @testset "sum(x; dims=$dims)" for dims in (:, 2, (1,3))
        # Forward
        test_frule(sum, rand(5); fkwargs=(;dims=dims))
        test_frule(sum, rand(ComplexF64, 2,3,4); fkwargs=(;dims=dims))

        # Reverse
        test_rrule(sum, rand(5); fkwargs=(;dims=dims))
        test_rrule(sum, rand(ComplexF64, 2,3,4); fkwargs=(;dims=dims))

        # Structured matrices
        test_rrule(sum, rand(5)'; fkwargs=(;dims=dims))
        y, back = rrule(sum, UpperTriangular(rand(5,5)); dims=dims)
        unthunk(back(y*(1+im))[2]) isa UpperTriangular{Float64}
        @test_skip test_rrule(sum, UpperTriangular(rand(5,5)) ⊢ randn(5,5); fkwargs=(;dims=dims), check_inferred=false) # Problem: in add!!  Evaluated: isapprox

        # Boolean -- via @non_differentiable
        test_rrule(sum, randn(5) .> 0; fkwargs=(;dims=dims))

        # Function allowing for 2nd derivatives
        for x in (rand(5), rand(2,3,4))
            dy = maximum(x; dims=dims)
            test_frule(ChainRules._unsum, x, dy, dims)
            test_rrule(ChainRules._unsum, x, dy, dims)
        end

        # Arrays of arrays
        for x  in ([rand(ComplexF64, 3) for _ in 1:4], [rand(3) for _ in 1:2, _ in 1:3, _ in 1:4])
            test_rrule(sum, x; fkwargs=(;dims=dims), check_inferred=false)

            dy = sum(x; dims=dims)
            ddy = rrule(ChainRules._unsum, x, dy, dims)[2](x)[3]
            @test size(ddy) == size(dy)
        end
    end

    @testset "sum abs2" begin
        sizes = (3, 4, 7)
        @testset "dims = $dims" for dims in (:, 1)
            @testset "Array{$N, $T}" for N in eachindex(sizes), T in (Float64, ComplexF64)
                x = randn(T, sizes[1:N]...)
                test_frule(sum, abs2, x; fkwargs=(;dims=dims))
                test_rrule(sum, abs2, x; fkwargs=(;dims=dims))
            end

            # Boolean -- via @non_differentiable, test that this isn't ambiguous
            test_rrule(sum, abs2, randn(5) .> 0; fkwargs=(;dims=dims))
        end
    end  # sum abs2

    @testset "sum(f, xs)" begin
        # This calls back into AD
        test_rrule(sum, abs, [-4.0, 2.0, 2.0])
        test_rrule(sum, cbrt, randn(5))
        test_rrule(sum, Multiplier(2.0), [2.0, 4.0, 8.0])

        # Complex numbers
        test_rrule(sum, sqrt, rand(ComplexF64, 5))
        test_rrule(sum, abs, rand(ComplexF64, 3, 4))  # complex -> real

        # inference fails for array of arrays
        test_rrule(sum, sum, [[2.0, 4.0], [4.0,1.9]]; check_inferred=false)

        # dims kwarg
        test_rrule(sum, abs, [-2.0 4.0; 5.0 1.9]; fkwargs=(;dims=1))
        test_rrule(sum, abs, [-2.0 4.0; 5.0 1.9]; fkwargs=(;dims=2))
        test_rrule(sum, sqrt, rand(ComplexF64, 3, 4); fkwargs=(;dims=(1,)))

        test_rrule(sum, abs, @SVector[1.0, -3.0])

        # covectors
        x = [-4.0 2.0; 2.0 -1.0]
        test_rrule(sum, inv, x[1, :]')
        test_rrule(sum, inv, x[1:1, :]')
        test_rrule(sum, inv, transpose(view(x, 1, :)))

        # Make sure we preserve type for StaticArrays
        _, pb = rrule(CFG, sum, abs, @SVector[1.0, -3.0])
        @test pb(1.0) isa Tuple{NoTangent, NoTangent, SVector{2, Float64}}

        # make sure we preserve type for Diagonal
        _, pb = rrule(CFG, sum, abs, Diagonal([1.0, -3.0]))
        @test pb(1.0)[3] isa Diagonal

        # Boolean -- via @non_differentiable, test that this isn't ambiguous
        test_rrule(sum, sqrt, randn(5) .> 0)
        test_rrule(sum, sqrt, randn(5,5) .> 0; fkwargs=(;dims=1))
        # ... and Bool produced by function
        @test_skip test_rrule(sum, iszero, randn(5))  # DimensionMismatch("second dimension of A, 1, does not match length of x, 0")

        # Functions that return a Vector
        # see https://github.com/FluxML/Zygote.jl/issues/1074
        test_rrule(sum, make_two_vec, [1.0, 3.0, 5.0, 7.0])
        test_rrule(sum, make_two_vec, [1.0 2.0; 3.0 4.0])
        test_rrule(sum, make_two_vec, [1.0 2.0; 3.0 4.0]; fkwargs=(;dims=2))
        test_rrule(sum, make_two_vec, [1.0 2.0; 3.0 4.0]; fkwargs=(;dims=1))
        test_rrule(sum, make_two_vec, [1.0 2.0; 3.0 4.0]; fkwargs=(;dims=(3, 4)))
    end

    # https://github.com/JuliaDiff/ChainRules.jl/issues/522
    @testset "sum(xs, weights) (#522)" begin
        xs = rand(5)
        weights = rand(5)

        @test rrule(SumRuleConfig(), Base.sum, xs, weights) isa Nothing
    end

    @testset "prod" begin
        @testset "Array{$T}" for T in [Float64, ComplexF64]
            @testset "size = $sz, dims = $dims" for (sz, dims) in [
                ((12,), :), ((12,), 1),
                ((3,4), 1), ((3,4), 2), ((3,4), :), ((3,4), [1,2]),
                ((3,4,1), 1), ((3,2,2), 3), ((3,2,2), 2:3),
                ]
                x = randn(T, sz)
                test_rrule(prod, x; fkwargs=(dims=dims,), check_inferred=true)
                x[1] = 0
                test_rrule(prod, x; fkwargs=(dims=dims,), check_inferred=true)
                x[5] = 0
                test_rrule(prod, x; fkwargs=(dims=dims,), check_inferred=true)
                x[3] = x[7] = 0  # two zeros along some slice, for any dims
                test_rrule(prod, x; fkwargs=(dims=dims,), check_inferred=true)

                if ndims(x) == 3
                    xp = PermutedDimsArray(x, (3,2,1))  # not a StridedArray
                    test_rrule(prod, xp; fkwargs=(dims=dims,), check_inferred=true)
                end
            end

            @testset "structured wrappers" begin
                # Adjoint -- like PermutedDimsArray this may actually be used
                xa = adjoint(rand(T,4,4))
                test_rrule(prod, xa)
                test_rrule(prod, xa, fkwargs=(dims=2,))
                # use Adjoint for tangent
                test_rrule(prod, xa ⊢ rand(T,4,4)')
                test_rrule(prod, xa ⊢ rand(T,4,4)', fkwargs=(dims=2,))
                @test unthunk(rrule(prod, adjoint(rand(T,3,3)))[2](1.0)[2]) isa Matrix
                @test unthunk(rrule(prod, adjoint(rand(T,3,3)), dims=1)[2](ones(1,3))[2]) isa Matrix

                # Diagonal -- a stupid thing to do, product of zeros! Shouldn't be an error though:
                @test iszero(unthunk(rrule(prod, Diagonal(rand(T,3)))[2](1.0)[2]))
                @test iszero(unthunk(rrule(prod, Diagonal(rand(T,3)), dims=1)[2](ones(1,3))[2]))
                # does a division for the complex case, so is not necessarily exact
                @test isapprox(
                    unthunk(rrule(prod, Diagonal(rand(T,1)))[2](1.0)[2]), # 1x1 sparse matrix
                    hcat(1);
                    rtol=T <: Complex ? 2eps() : 0.0,
                )
                @test unthunk(rrule(prod, Diagonal(ones(T,2)), dims=1)[2](ones(1,2))[2]) == Diagonal([0.0 1; 1 0])

                # Triangular -- almost equally stupud
                @test iszero(unthunk(rrule(prod, UpperTriangular(rand(T,3,3)))[2](1.0)[2]))
                @test unthunk(rrule(prod, UpperTriangular(ones(T,2,2)))[2](1.0)[2]) == UpperTriangular([0.0 0; 1 0])

                # Symmetric -- at least this doesn't have zeros, still an unlikely combination
                xs = Symmetric(rand(T,4,4))
                @test unthunk(rrule(prod, Symmetric(T[1 2; -333 4]))[2](1.0)[2]) == [16 8; 8 4]
                # TODO debug why these fail  https://github.com/JuliaDiff/ChainRules.jl/issues/475
                @test_skip test_rrule(prod, xs)
                @test_skip test_rrule(prod, xs, fkwargs=(dims=2,))
            end
        end
        @testset "Array{Float32}, no zero entries" begin
            v = [1f-5, 1f-10, 1f-15, 1f-20]
            @test prod(v) == 0
            @test unthunk(rrule(prod, v)[2](1f0)[2]) == zeros(4)
            test_rrule(prod, v)
        end
    end  # prod

    @testset "foldl(f, ::Array)" begin
        # Simple
        y1, b1 = rrule(CFG, foldl, *, [1, 2, 3]; init=1)
        @test y1 == 6
        b1(7) == (NoTangent(), NoTangent(), [42, 21, 14])

        y2, b2 = rrule(CFG, foldl, *, [1 2; 0 4])  # without init, needs vcat
        @test y2 == 0
        b2(8) == (NoTangent(), NoTangent(), [0 0; 64 0])  # matrix, needs reshape

        # Test execution order
        c5 = Counter()
        y5, b5 = rrule(CFG, foldl, c5, [5, 7, 11])
        @test c5 == Counter(2)
        @test y5 == ((5 + 7)*1 + 11)*2 == foldl(Counter(), [5, 7, 11])
        @test b5(1) == (NoTangent(), NoTangent(), [12*32, 12*42, 22])
        @test c5 == Counter(42)

        c6 = Counter()
        y6, b6 = rrule(CFG, foldl, c6, [5, 7, 11], init=3)
        @test c6 == Counter(3)
        @test y6 == (((3 + 5)*1 + 7)*2 + 11)*3 == foldl(Counter(), [5, 7, 11], init=3)
        @test b6(1) == (NoTangent(), NoTangent(), [63*33*13, 43*13, 23])
        @test c6 == Counter(63)

        # Test gradient of function
        y7, b7 = rrule(CFG, foldl, Multiplier(3), [5, 7, 11])
        @test y7 == foldl((x,y)->x*y*3, [5, 7, 11])
        @test b7(1) == (NoTangent(), Tangent{Multiplier{Int}}(x = 2310,), [693, 495, 315])

        y8, b8 = rrule(CFG, foldl, Multiplier(13), [5, 7, 11], init=3)
        @test y8 == 2_537_535 == foldl((x,y)->x*y*13, [5, 7, 11], init=3)
        @test b8(1) == (NoTangent(), Tangent{Multiplier{Int}}(x = 585585,), [507507, 362505, 230685])
        # To find these numbers:
        # ForwardDiff.derivative(z -> foldl((x,y)->x*y*z, [5,7,11], init=3), 13)
        # ForwardDiff.gradient(z -> foldl((x,y)->x*y*13, z, init=3), [5,7,11]) |> string

        # Finite differencing
        test_rrule(foldl, /, 1 .+ rand(3,4))
        test_rrule(foldl, *, rand(ComplexF64,3,4); fkwargs=(; init=rand(ComplexF64)))
        test_rrule(foldl, +, rand(ComplexF64,7); fkwargs=(; init=rand(ComplexF64)))
        test_rrule(foldl, max, rand(3); fkwargs=(; init=999))
    end
    VERSION >= v"1.5" && @testset "foldl(f, ::Tuple)" begin
        y1, b1 = rrule(CFG, foldl, *, (1,2,3); init=1)
        @test y1 == 6
        b1(7) == (NoTangent(), NoTangent(), Tangent{NTuple{3,Int}}(42, 21, 14))

        y2, b2 = rrule(CFG, foldl, *, (1, 2, 0, 4))
        @test y2 == 0
        b2(8) == (NoTangent(), NoTangent(), Tangent{NTuple{4,Int}}(0, 0, 64, 0))

        # Finite differencing
        test_rrule(foldl, /, Tuple(1 .+ rand(5)))
        test_rrule(foldl, *, Tuple(rand(ComplexF64, 5)))
    end
end

@testset "Accumulations" begin
    @testset "cumprod" begin
        v = round.(10 .* randn(9), sigdigits=3)
        test_rrule(cumprod, v)
        v[3] = 0
        test_rrule(cumprod, v)
        v[6] = 0
        test_rrule(cumprod, v)

        @testset "higher dimensions, dims=$dims" for dims in (1,2,3)
            m = round.(10 .* randn(4,5), sigdigits=3)
            test_rrule(cumprod, m; fkwargs=(;dims=dims), atol=0.1)
            m[2, 2] = 0
            m[2, 4] = 0
            test_rrule(cumprod, m; fkwargs=(;dims=dims))

            t = round.(10 .* randn(3,3,3), sigdigits=3)
            test_rrule(cumprod, t; fkwargs=(;dims=dims))
            t[2, 2, 2] = 0
            t[2, 3, 3] = 0
            test_rrule(cumprod, t; fkwargs=(;dims=dims))
        end

        @testset "types" begin
            back = rrule(cumprod, [1, 2, 3])[2]  # rule allows integer input, but test_rrule does not
            @test unthunk(back(fill(0.5, 3))[2]) == [9/2, 2, 1]

            back = rrule(cumprod, PermutedDimsArray([1 2; 3 4], (2,1)); dims=1)[2]
            @test unthunk(back(ones(Float32, 2,2))[2]) == [3 5; 1 3]

            @test_throws Exception cumprod(Symmetric([1 2; 3 4]), dims=1) # forward pass fails, so can't test gradient

            back = rrule(cumprod, Diagonal([1, 2]); dims=1)[2]
            @test unthunk(back(fill(0.5, 2, 2))[2]) ≈ [1/2 0; 0 0]  # ProjectTo'd to Diagonal now
        end
    end  # cumprod

    @testset "accumulate(f, ::Array)" begin
        # Simple
        y1, b1 = rrule(CFG, accumulate, *, [1, 2, 3, 4]; init=1)
        @test y1 == [1, 2, 6, 24]
        @test b1([1, 1, 1, 1]) == (NoTangent(), NoTangent(), [33, 16, 10, 6])

        if VERSION >= v"1.5"
            y2, b2 = rrule(CFG, accumulate, /, [1 2; 3 4])
            @test y2 ≈ accumulate(/, [1 2; 3 4])
            @test b2(ones(2, 2))[3] ≈ [1.5416666 -0.104166664; -0.18055555 -0.010416667]  atol=1e-6
        end

        # Test execution order
        c3 = Counter()
        y3, b3 = rrule(CFG, accumulate, c3, [5, 7, 11]; init=3)
        @test c3 == Counter(3)
        @test y3 == [8, 30, 123] == accumulate(Counter(), [5, 7, 11]; init=3)
        @test b3([1, 1, 1]) == (NoTangent(), NoTangent(), [29169, 602, 23]) # the 23 is clear!

        c4 = Counter()
        y4, b4 = rrule(CFG, accumulate, c4, [5, 7, 11])
        @test c4 == Counter(2)
        @test y4 == [5, (5+7)*1, ((5+7)*1 + 11)*2] == accumulate(Counter(), [5, 7, 11])
        @test b4([1, 1, 1]) == (NoTangent(), NoTangent(), [417, 42*(1 + 12), 22])

        # Test gradient of function
        y7, b7 = rrule(CFG, accumulate, Multiplier(3), [5, 7, 11])
        @test y7 == accumulate((x,y)->x*y*3, [5, 7, 11])
        @test b7([1, 1, 1]) == (NoTangent(), Tangent{Multiplier{Int}}(x = 2345,), [715, 510, 315])

        y8, b8 = rrule(CFG, accumulate, Multiplier(13), [5, 7, 11], init=3)
        @test y8 == [195, 17745, 2537535] == accumulate((x,y)->x*y*13, [5, 7, 11], init=3)
        @test b8([1, 1, 1]) == (NoTangent(), Tangent{Multiplier{Int}}(x = 588330,), [511095, 365040, 230685])
        # To find these numbers:
        # ForwardDiff.derivative(z -> sum(accumulate((x,y)->x*y*z, [5,7,11], init=3)), 13)
        # ForwardDiff.gradient(z -> sum(accumulate((x,y)->x*y*13, z, init=3)), [5,7,11]) |> string

        # Finite differencing
        test_rrule(accumulate, *, randn(5); fkwargs=(; init=rand()))
        if VERSION >= v"1.5"
            test_rrule(accumulate, /, 1 .+ rand(3, 4))
            test_rrule(accumulate, ^, 1 .+ rand(2, 3); fkwargs=(; init=rand()))
        end
    end
    VERSION >= v"1.5" && @testset "accumulate(f, ::Tuple)" begin
        # Simple
        y1, b1 = rrule(CFG, accumulate, *, (1, 2, 3, 4); init=1)
        @test y1 == (1, 2, 6, 24)
        @test b1((1, 1, 1, 1)) == (NoTangent(), NoTangent(), Tangent{NTuple{4,Int}}(33, 16, 10, 6))

        # Finite differencing
        test_rrule(accumulate, *, Tuple(randn(5)); fkwargs=(; init=rand()))
        test_rrule(accumulate, /, Tuple(1 .+ rand(5)); check_inferred=false)
    end
end


text:
using Gaugefields
using Documenter

DocMeta.setdocmeta!(Gaugefields, :DocTestSetup, :(using Gaugefields); recursive=true)

makedocs(;
    modules=[Gaugefields],
    authors="Akio Tomiya, Yuki Nagai <cometscome@gmail.com> and contributors",
    repo="https://github.com/akio-tomiya/Gaugefields.jl/blob/{commit}{path}#{line}",
    sitename="Gaugefields.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://github.com/akio-tomiya/Gaugefields.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        #"File loading" => "fileloading.md",
        #"Heatbath updates" => "heatbath.md",
        #"Gradientflow" => "gradientflow.md",
        "Utilities" => "utilities.md",
        "Applications" => "applications.md",
        "Useful functions" => "usefulfunctions.md",
        "How to implement new gauge fields" => "howtoimplement.md",
        "Measurements" => "measurements.md",
        "Parallel computation" => "mpi.md",
        #"Derivatives" => "derivatives.md",
        #"Hybrid Monte Carlo" => "hmc.md",
        #"Smearing" => "smearing.md",
    ],
)

deploydocs(;
    repo="github.com/akio-tomiya/Gaugefields.jl",
    devbranch="master",
)


text:
# Benchmark script for comparing the convergence rate of GD, BiasSEGA, and SEGA
# for minimizing a multivariate quadratic on the unit ball

using Random, LinearAlgebra

function main(n=10, γ=1/(2n), niterations=30)

    Random.seed!(123)
    Q = Symmetric(randn(n, n))
    c = randn(n)
    f = (x) -> x'*Q*x ./ 2 .+ c'*x
    ∇ = (x) -> Q*x .+ c

    xgd = zeros(n, niterations)
    xbias = zeros(n, niterations)
    xsega = zeros(n, niterations)

    bsega = BiasSEGA(n)
    sega = SEGA(n)

    s = zeros(n)
    for i in 2:niterations
        # GD
        xgd[:, i] .= xgd[:, i-1] .- γ.*∇(xgd[:, i-1])
        xgd[:, i] ./= norm(xgd[:, i])

        # BiasSEGA
        s .= randn(n)
        project!(bsega, dot(s, ∇(xbias[:, i-1])), s)
        xbias[:, i] .= xbias[:, i-1] .- γ.*gradient(bsega)
        xbias[:, i] ./= norm(xbias[:, i])

        # SEGA
        project!(sega, 1/n, dot(s, ∇(xsega[:, i-1])), s)
        unbias!(sega)
        xsega[:, i] .= xsega[:, i-1] .- γ.*gradient(sega)
        xsega[:, i] ./= norm(xsega[:, i])
    end

    println("GD \t BiasSEGA \t SEGA")
    for i in 2:niterations
        println("$(f(xgd[:, i])) \t $(f(xbias[:, i])) \t $(f(xsega[:, i]))")
    end
end


text:
include("../../testcases.jl")
using DFTK

atoms = [ElementPsp(:Si, psp=load_psp("hgh/pbe/si-q4")) => silicon.positions]
model = model_DFT(silicon.lattice, atoms, [:gga_c_pbe, :gga_x_pbe])

abinitpseudos = [joinpath(@__DIR__, "Si-q4-pbe.abinit.hgh")]
DFTK.run_abinit_scf(model, @__DIR__; abinitpseudos=abinitpseudos,
                    Ecut=25, kgrid=[4, 4, 4], n_bands=6, tol=1e-10,
                    iscf=3) # Use Anderson mixing instead of minimization


text:
# ---
# title: Histograms
# cover: assets/histograms.png
# description: Computing 1- and 2-dimensional histograms
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
using AlgebraOfGraphics: density
set_aog_theme!() #src

df = (x=rand(0:99, 1000),)
plt = data(df) * mapping(:x) * histogram(bins=20)
fg = draw(plt)

#

df = (x=randn(1000), c=rand(["a", "b"], 1000))
plt = data(df) * mapping(:x, color=:c, stack=:c) * histogram(bins=20)
fg = draw(plt)

#

df = (x=rand(1000), y=randn(1000))
plt = data(df) * mapping(:x, :y) * histogram(bins=20)
draw(plt)

# save cover image #src
mkpath("assets") #src
save("assets/histograms.png", fg) #src


text:
using LinearAlgebra, MATLAB, ForwardDiff, SparseArrays, Printf
using Convex, Hypatia
const FD = ForwardDiff
include("/Users/kevintracy/.julia/dev/CPEG/src/qp_solver.jl")

function dynamics(x,u)
    r = x[1:3]
    v = x[4:6]
    return [v;u[1:3] + [0;0;-9.81]]
end
function rk2(x,u)
    Δt = u[4]
    k1 = Δt*dynamics(x,u)
    k2 = Δt*dynamics(x + k1/2, u)
    return x + k2
end
function rollout(x0,U)
    N = length(U)+1
    X = [zeros(length(x0)) for i = 1:N]
    X[1] = x0
    for i = 1:(N-1)
        X[i+1] = rk2(X[i],U[i])
    end
    return X
end
function get_jacobians(X,U)
    N = length(X)
    A = [FD.jacobian(_x -> rk2(_x,U[i]),X[i]) for i = 1:N-1]
    B = [FD.jacobian(_u -> rk2(X[i],_u),U[i]) for i = 1:N-1]
    return A,B
end

function mpc(X,U)
    A,B = get_jacobians(X,U)
    nx = length(X[1])
    nu = length(U[1])
    N = length(X)
    alt_target = 10.0
    dt_target = 0.5

    δX = Variable(nx,N)
    δU = Variable(nu,N-1)

    J = 0

    γ = 1e3
    for i = 1:N-1
        J += γ*quadform(U[i][1:3] + δU[1:3,i],Matrix(I(3)))
        J += square(U[i][4] + δU[4,i] - dt_target)
    end

    p= minimize(J)

    for i = 1:N-1
        p.constraints += δU[4,i] >= 0.3
        p.constraints += δU[4,i] <= 0.7
    end
    p.constraints += X[N][3] + δX[3,N] == alt_target
    p.constraints += δX[:,1] == zeros(6)
    for i = 1:N-1
        # @info "here?"
        p.constraints += δX[:,i+1] == A[i]*δX[:,i] + B[i]*δU[:,i]
        # @info "nope"
    end

    solve!(p, Hypatia.Optimizer();silent_solver = true)

    δU = evaluate(δU)
    # @show norm(δU)
    U = [U[i] + vec(δU[:,i]) for i = 1:length(U)]

    return U, norm(δU)
end

function mat_from_vec(X)
    Xm = zeros(length(X[1]),length(X))
    for i = 1:length(X)
        Xm[:,i] = X[i]
    end
    return Xm
end
function tt()

    x0 = [0;0;100;0;0;-10.0]
    # u0 = [0;0;9.7]
    u1 = [0;0;9.81;0.5]

    N = 20

    U = [copy(u1) for i = 1:N-1]

    X = rollout(x0,U)

    xm = mat_from_vec(X)

    for i = 1:5

        U, ndu = mpc(X,U)

        if ndu < 1e-1
            break
        end
        X = rollout(x0,U)
        @show X[end][3]

    end

    mat"
    figure
    hold on
    %plot3($xm(1,:),$xm(2,:),$xm(3,:))
    plot($xm(1:3,:)')
    hold off
    "

    return nothing
end

tt()


text:
include("pot_id.jl")
include("pot.jl")
include("message_handling.jl")
include("scheduling.jl")