text:
"""
Type: NLPAtX

Methods: update!, reinit!

NLPAtX contains the information concerning a nonlinear optimization model at
the iterate x.

min_{x ∈ ℜⁿ} f(x) subject to lcon <= c(x) <= ucon, lvar <= x <= uvar.

Tracked data include:
 - x             : the current iterate
 - fx [opt]      : function evaluation at x
 - gx [opt]      : gradient evaluation at x
 - Hx [opt]      : hessian evaluation at x

 - mu [opt]      : Lagrange multiplier of the bounds constraints

 - cx [opt]      : evaluation of the constraint function at x
 - Jx [opt]      : jacobian matrix of the constraint function at x
 - lambda        : Lagrange multiplier of the constraints

 - d [opt]       : search direction
 - res [opt]     : residual

 - current_time  : time
 - current_score : score
 (import the type NLPModels.Counters)

Constructors:
 `NLPAtX(:: T, :: T, :: S; fx :: eltype(T) = _init_field(eltype(T)), gx :: T = _init_field(T), Hx :: Matrix{eltype(T)} = _init_field(Matrix{eltype(T)}), mu :: T = _init_field(T), cx :: T = _init_field(T), Jx :: Matrix{eltype(T)} = _init_field(Matrix{eltype(T)}), d :: T = _init_field(T), res :: T = _init_field(T), current_time :: Float64 = NaN) where {S, T <: AbstractVector}`

 `NLPAtX(:: T; fx :: eltype(T) = _init_field(eltype(T)), gx :: T = _init_field(T), Hx :: Matrix{eltype(T)} = _init_field(Matrix{eltype(T)}), mu :: T = _init_field(T), current_time :: Float64 = NaN, current_score :: Union{T,eltype(T)} = _init_field(eltype(T))) where {T <: AbstractVector}`

 `NLPAtX(:: T, :: T; fx :: eltype(T) = _init_field(eltype(T)), gx :: T = _init_field(T), Hx :: Matrix{eltype(T)} = _init_field(Matrix{eltype(T)}), mu :: T = _init_field(T), cx :: T = _init_field(T), Jx :: Matrix{eltype(T)} = _init_field(Matrix{eltype(T)}), d :: T = _init_field(T), res :: T = _init_field(T), current_time :: Float64  = NaN, current_score :: Union{T,eltype(T)} = _init_field(eltype(T))) where T <: AbstractVector`

Note:
 - By default, unknown entries are set using `_init_field`.
 - By default the type of `current_score` is `eltype(x)` and cannot be changed once the State is created.
    To have a vectorized `current_score` of length n, try something like `GenericState(x, Array{eltype(x),1}(undef, n))`.
 - All these information (except for `x` and `lambda`) are optionnal and need to be update when
    required. The update is done through the `update!` function.
 - `x` and `lambda` are mandatory entries. If no constraints `lambda = []`.
 - The constructor check the size of the entries.

See also: `GenericState`, `update!`, `update_and_start!`, `update_and_stop!`, `reinit!`
"""
mutable struct 	NLPAtX{S,
                       T  <: AbstractVector,
                       HT <: AbstractMatrix,
                       JT <: AbstractMatrix}  <: AbstractState{S, T}

  #Unconstrained State
  x            :: T     # current point
  fx           :: eltype(T) # objective function
  gx           :: T  # gradient size: x
  Hx           :: HT  # hessian size: |x| x |x|

  #Bounds State
  mu           :: T # Lagrange multipliers with bounds size of |x|

  #Constrained State
  cx           :: T # vector of constraints lc <= c(x) <= uc
  Jx           :: JT  # jacobian matrix, size: |lambda| x |x|
  lambda       :: T    # Lagrange multipliers

  d            :: T #search direction
  res          :: T #residual

  #Resources State
  current_time   :: Float64
  current_score  :: S

  function NLPAtX(x            :: T,
                 lambda        :: T,
                 current_score :: S;
                 fx            :: eltype(T) = _init_field(eltype(T)),
                 gx            :: T = _init_field(T),
                 Hx            :: AbstractMatrix = _init_field(Matrix{eltype(T)}),
                 mu            :: T = _init_field(T),
                 cx            :: T = _init_field(T),
                 Jx            :: AbstractMatrix = _init_field(Matrix{eltype(T)}),
                 d             :: T = _init_field(T),
                 res           :: T = _init_field(T),
                 current_time  :: Float64 = NaN,
                 ) where {S, T <: AbstractVector}

    _size_check(x, lambda, fx, gx, Hx, mu, cx, Jx)

    return new{S, T, typeof(Hx), typeof(Jx)}(
      x,
      fx,
      gx,
      Hx,
      mu,
      cx,
      Jx,
      lambda,
      d,
      res,
      current_time,
      current_score,
    )
  end
end

function NLPAtX(x             :: T,
                lambda        :: T;
                fx            :: eltype(T) = _init_field(eltype(T)),
                gx            :: T = _init_field(T),
                Hx            :: AbstractMatrix = _init_field(Matrix{eltype(T)}),
                mu            :: T = _init_field(T),
                cx            :: T = _init_field(T),
                Jx            :: AbstractMatrix = _init_field(Matrix{eltype(T)}),
                d             :: T = _init_field(T),
                res           :: T = _init_field(T),
                current_time  :: Float64  = NaN,
                current_score :: Union{T,eltype(T)} = _init_field(eltype(T)),
                ) where T <: AbstractVector

  _size_check(x, lambda, fx, gx, Hx, mu, cx, Jx)

  return NLPAtX(
    x, lambda, current_score,
    fx = fx, gx = gx, Hx = Hx,
    mu = mu, cx = cx, Jx = Jx,
    d = d, res = res,
    current_time = current_time,
  )
end

function NLPAtX(x             :: T;
                fx            :: eltype(T) = _init_field(eltype(T)),
                gx            :: T = _init_field(T),
                Hx            :: AbstractMatrix = _init_field(Matrix{eltype(T)}),
                mu            :: T = _init_field(T),
                d             :: T = _init_field(T),
                res           :: T = _init_field(T),
                current_time  :: Float64 = NaN,
                current_score :: Union{T,eltype(T)} = _init_field(eltype(T)),
                ) where {T <: AbstractVector}

  _size_check(x, zeros(eltype(T),0), fx, gx, Hx, mu,
                _init_field(T), _init_field(Matrix{eltype(T)}))

  return NLPAtX(x, zeros(eltype(T),0), current_score,
                fx = fx, gx = gx,
                Hx = Hx, mu = mu,
                d = d, res = res,
                current_time = current_time)
end

"""
reinit!: function that set all the entries at void except the mandatory x

`reinit!(:: NLPAtX, x :: AbstractVector, l :: AbstractVector; kwargs...)`

`reinit!(:: NLPAtX; kwargs...)`

Note: if `x` or `lambda` are given as keyword arguments they will be
prioritized over the existing `x`, `lambda` and the default `Counters`.
"""
function reinit!(stateatx :: NLPAtX{S, T, MT},
                 x        :: T,
                 l        :: T;
                 kwargs...) where {S, T, MT}

  for k ∈ fieldnames(NLPAtX)
    if k ∉ [:x,:lambda]
      setfield!(stateatx, k, _init_field(typeof(getfield(stateatx, k))))
    end
  end

  setfield!(stateatx, :x, x)
  setfield!(stateatx, :lambda, l)

  if length(kwargs)==0
    return stateatx #save the update! call if no other kwargs than x
  end

  return update!(stateatx; kwargs...)
end

function reinit!(stateatx :: NLPAtX{S, T, MT},
                 x        :: T;
                 kwargs...) where {S, T, MT}

  for k ∈ fieldnames(NLPAtX)
    if k ∉ [:x,:lambda]
      setfield!(stateatx, k, _init_field(typeof(getfield(stateatx, k))))
    end
  end

  setfield!(stateatx, :x, x)

  if length(kwargs)==0
    return stateatx #save the update! call if no other kwargs than x
  end

  return update!(stateatx; kwargs...)
end

function reinit!(stateatx :: NLPAtX; kwargs...)

  for k ∈ fieldnames(NLPAtX)
    if k ∉ [:x,:lambda]
      setfield!(stateatx, k, _init_field(typeof(getfield(stateatx, k))))
    end
  end

  return update!(stateatx; kwargs...)
end

"""
_size_check!: check the size of the entries in the State

`_size_check(x, lambda, fx, gx, Hx, mu, cx, Jx)`
"""
function _size_check(x, lambda, fx, gx, Hx, mu, cx, Jx)

  if length(gx) != 0 &&  length(gx) != length(x)
    throw(error("Wrong size of gx in the NLPAtX."))
  end
  if size(Hx) != (0,0) && size(Hx) != (length(x), length(x))
    throw(error("Wrong size of Hx in the NLPAtX."))
  end
  if length(mu) != 0 && length(mu) != length(x)
    throw(error("Wrong size of mu in the NLPAtX."))
  end

  if lambda != zeros(0)
    if length(cx) != 0 && length(cx) != length(lambda)
      throw(error("Wrong size of cx in the NLPAtX."))
    end
    if size(Jx) != (0,0) && size(Jx) != (length(lambda), length(x))
      throw(error("Wrong size of Jx in the NLPAtX."))
    end
  end
end


text:
#
using PyPlot, DelimitedFiles, DSP, Statistics, FFTW

if !isdefined(Main, :Ganymede)
   include("/Users/hyzhou/Documents/Computer/Julia/BATSRUS/VisAnaJulia/space/Ganymede.jl")
   using .Ganymede: Rg, upstream_value
end

function smooth(x::Vector, n=100)
   nx = length(x)
   x̄ = zeros(eltype(x),length(x))

   # left points
   for i = 1:n
      x̄[i] = mean(x[1:(i+n)])
   end

   # middle points
   for i = (n+1):(nx-n)
      x̄[i] = mean(x[(i-n):(i+n)])
   end

   # right points
   for i = (nx-n+1):nx
      x̄[i] = mean(x[(i-n):end])
   end

   return x̄
end

# Background
# G8 B= -10 -6 -86 [nT], U=140 0 0 [km/s]
Uxbk = upstream_value[2]
Bybk = upstream_value[6]
Bzbk = upstream_value[7]
TiltedAngle = atan(Bybk/Bzbk)

Δy = 5.64 # [Rg]

# Take clock angle in yz plane into consideration
Potential_bk = Δy * cos(TiltedAngle) * Rg *
   Uxbk*1e3 * abs(Bzbk*cos(TiltedAngle)+Bybk*sin(TiltedAngle))*1e-9 *
   1e-3 # [kV]

dir = "/Users/hyzhou/Documents/Computer/ParaView/data/"
dirHallBoundary = "dataBoundary_Hall/"
dirPICBoundary = "dataBoundary_PIC/"
dirHallSurface = "dataSurf_Hall/"
dirPICSurface = "dataSurf_PIC/"

data_hall = readdlm(dir*"CPCP_hall.txt", comments=true)
data_pic  = readdlm(dir*"CPCP_pic.txt", comments=true)

data_hall_b = readdlm(dirHallBoundary*"CPCP_test.txt", comments=true)
data_pic_b  = readdlm(dirPICBoundary*"CPCP_test.txt", comments=true)

data_hall_s = readdlm(dirHallSurface*"CPCP_test.txt", comments=true)
data_pic_s  = readdlm(dirPICSurface*"CPCP_test.txt", comments=true)

# Check missing frames
#missingtime = setdiff(0:1200,data_hall[:,1])
# hall: [355, 457, 505, 568]

#missingtime = setdiff(300:1500,data_pic[:,1])
# pic: [567, 591, 610, 680]

Fs = 1                 # Sampling frequency
T = 1/Fs               # Sampling period
L = size(data_hall_b)[1] # Length of signal
t = (0:L-1)*T          # Time vector

#missingtime = setdiff(0:1200,data_hall[:,1]) # Check missing frames

e_hall = data_hall[:,2] ./ Potential_bk
e_pic  = data_pic[:,2] ./ Potential_bk
e_hall_mean = mean(e_hall)
e_hall_std = std(e_hall)
e_pic_mean = mean(e_pic)
e_pic_std = std(e_pic)

e_hall_b = data_hall_b[:,2] ./ Potential_bk
e_pic_b  = data_pic_b[:,2] ./ Potential_bk
e_hall_mean_b = mean(e_hall_b)
e_hall_std_b = std(e_hall_b)
e_pic_mean_b = mean(e_pic_b)
e_pic_std_b = std(e_pic_b)

e_hall_s = data_hall_s[:,2] ./ Potential_bk
e_pic_s  = data_pic_s[:,2] ./ Potential_bk

#=
# Strange points: hall -> 13s, 14s, 415s, PIC -> 1255s, 1256s
# Remove strange points
for i = 2:length(e_hall_s)
   if e_hall_s[i] > 1
      e_hall_s[i] = 1.
   elseif e_hall_s[i] < 0
      e_hall_s[i] = 0.
   end
   if e_hall_s[i] - e_hall_s[i-1] > 0.2
      e_hall_s[i] = 0.5 * (e_hall_s[i] + e_hall_s[i-1])
   end
end

for i = 2:length(e_pic_s)
   if e_pic_s[i] > 1 || e_pic_s[i] < 0
      e_pic_s[i] = 0.0
   end
   if e_pic_s[i] - e_pic_s[i-1] > 0.2
      e_pic_s[i] = 0.5 * (e_pic_s[i] + e_pic_s[i-1])
   end
end
=#

e_hall_s = smooth(e_hall_s, 2)
e_pic_s = smooth(e_pic_s, 2)

e_hall_mean_s = mean(e_hall_s)
e_hall_std_s = std(e_hall_s)
e_pic_mean_s = mean(e_pic_s)
e_pic_std_s = std(e_pic_s)

fig, ax = plt.subplots(2,1,figsize=(10.0,5.5))
plt.rc("font", family = "serif", size = 13)

#ax[1].plot(data_hall[:,1], e_hall, label="Hall MHD Old")
ax[1].plot(data_hall_b[:,1], e_hall_b, label="Hall MHD")
ax[1].plot(data_hall_s[:,1], e_hall_s, alpha=0.8, label="Hall MHD S")

ax[1].legend()

#ax[2].plot(data_pic[:,1].+300, e_pic, label="MHD-EPIC Old")
ax[2].plot(data_pic_b[:,1], e_pic_b, label="MHD-EPIC")
ax[2].plot(data_pic_s[:,1], e_pic_s, alpha=0.8, label="MHD-EPIC S")

ax[2].legend()

Yhall = welch_pgram(data_hall_s[:,2], 180, 100; onesided=true, fs=Fs)
Ypic = welch_pgram(data_pic_s[:,2], 180, 100; onesided=true, fs=Fs)

Yhall_fft = fft(data_hall_s[:,2])
Ypic_fft = fft(data_pic_s[:,2])

P2_hall = @. abs(Yhall_fft/L)
P1_hall = P2_hall[1:floor(Int,L/2+1)]
P1_hall[2:end-1] .= 2 .* P1_hall[2:end-1]

f_hall = Fs*(0:(L/2))/L

P2_pic = @. abs(Ypic_fft/L)
P1_pic = P2_pic[1:floor(Int,L/2+1)]
P1_pic[2:end-1] .= 2 .* P1_pic[2:end-1]

f_pic = Fs*(0:(L/2))/L

fig, ax = plt.subplots(2,1,figsize=(10.0,5.5))
#plt.rc("font", family="serif", size=12)
ax[1].plot(1 ./ f_hall[2:end], P1_hall[2:end], label="Hall MHD")
ax[2].plot(1 ./ f_pic[2:end], P1_pic[2:end], label="MHD-EPIC")
ax[2].set_xlabel("Periods [s]")
#ax[2].set_ylabel("Power")
ax[1].grid("on")
ax[2].grid("on")
#plt.yscale("log")
#plt.xscale("log")
ax[2].legend()
ax[2].tick_params(which="both", direction="in")
ax[2].annotate("(b)", xy=(-0.03, 1.0), xycoords="axes fraction", fontsize=14)
tight_layout()

fig, ax = plt.subplots(2,1,figsize=(10.0,5.5))
#plt.rc("font", family="serif", size=12)
ax[1].plot(1 ./ Yhall.freq[2:end], Yhall.power[2:end], label="Hall MHD")
ax[2].plot(1 ./ Ypic.freq[2:end], Ypic.power[2:end], label="MHD-EPIC")
ax[2].set_xlabel("Periods [s]")
ax[2].set_ylabel("Power")
ax[1].grid("on")
ax[2].grid("on")
#plt.yscale("log")
#plt.xscale("log")
ax[2].legend()
ax[2].tick_params(which="both", direction="in")
ax[2].annotate("(b)", xy=(-0.03, 1.0), xycoords="axes fraction", fontsize=14)
tight_layout()

dirHallBoundary = "dataBoundary_Hall/"
dirPICBoundary = "dataBoundary_PIC/"

data_hall_b = readdlm(dirHallBoundary*"CPCP_test.txt", comments=true)
data_pic_b  = readdlm(dirPICBoundary*"CPCP_test.txt", comments=true)

data_test_hall = readdlm("dataSurf_Hall/CPCP_test.txt", comments=true)
data_test_pic = readdlm("dataSurf_PIC/CPCP_test.txt", comments=true)

figure()
plot(data_test_hall[:,1],data_test_hall[:,2], label="1:total of 2 and 3 ")
plot(data_test_hall[:,1],data_test_hall[:,3], label="2:time derivative of surface integral")
plot(data_test_hall[:,1],data_test_hall[:,4], label="3:middle line")
plot(data_hall_b[:,1], data_hall_b[:,2], label="4:upstream curve integral")
xlabel("time [s]")
ylabel("potential [kV]")
legend()
title("Hall MHD")
tight_layout()

figure()
plot(data_test_pic[:,1],data_test_pic[:,2], label="1:total of 2 and 3 ")
plot(data_test_pic[:,1],data_test_pic[:,3], label="2:time derivative of surface integral")
plot(data_test_pic[:,1],data_test_pic[:,4], label="3:middle line")
plot(data_pic_b[:,1], data_pic_b[:,2], label="4:upstream curve integral")
xlabel("time [s]")
ylabel("potential [kV]")
legend()
title("MHD-EPIC")
tight_layout()

ext: jl


text:
function load_documenter_docs(pkg)
  docs = docsdir(pkg)
  isempty(docs) && return []

  getfiles(docs)
end

function getfiles(path, files = Tuple{String, String}[])
  isdir(path) || return files
  for f in readdir(path)
    f = joinpath(path, f)
    if isfile(f) && split(f, '.')[end] == "md"
      push!(files, (f, read(f)))
    elseif isdir(f)
      getfiles(f, files)
    end
  end
  files
end

function searchfiles(needle, files::Vector{Tuple{String, String}})
  scores = Float64[]

  for (path, content) in files
    push!(scores, compare(TokenSet(Jaro()), needle, content))
  end
  p = sortperm(scores, rev=true)[1:min(20, length(scores))]
  scores[p], files[p]
end

searchfiles(needle, pkg::String) = searchfiles(load_documenter_docs(pkg))

function searchfiles(needle)
  files = Tuple{String, String}[]
  for pkg in readdir(Pkg.dir())
    isdir(Pkg.dir(pkg)) || continue
    append!(files, load_documenter_docs(pkg))
  end
  searchfiles(needle, files)
end


text:
"""
$(TYPEDEF)

A cone of squares on a Jordan algebra.
"""
abstract type ConeOfSquares{T <: Real} end

# a cache for a cone of squares implementation
abstract type CSqrCache{T <: Real} end

"""
$(TYPEDEF)

A univariate convex function defined on positive reals.
"""
abstract type SepSpectralFun end

"""
$(TYPEDEF)

Epigraph of perspective function of a convex separable spectral function `h`
over a cone of squares `Q` on a Jordan algebra with rank `d`.

    $(FUNCTIONNAME){Q, T}(h::Hypatia.Cones.SepSpectralFun, d::Int, use_dual::Bool = false)
"""
mutable struct EpiPerSepSpectral{Q <: ConeOfSquares, T <: Real} <: Cone{T}
    h::SepSpectralFun
    use_dual_barrier::Bool
    d::Int
    dim::Int

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_aux_updated::Bool
    inv_hess_aux_updated::Bool
    dder3_aux_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    w_view::SubArray{T, 1}
    cache::CSqrCache{T}

    function EpiPerSepSpectral{Q, T}(
        h::SepSpectralFun,
        d::Int; # dimension/rank parametrizing the cone of squares
        use_dual::Bool = false,
        ) where {T <: Real, Q <: ConeOfSquares{T}}
        @assert d >= 1
        cone = new{Q, T}()
        cone.h = h
        cone.use_dual_barrier = use_dual
        cone.d = d
        cone.dim = 2 + vector_dim(Q, d)
        return cone
    end
end

reset_data(cone::EpiPerSepSpectral) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.hess_aux_updated =
    cone.inv_hess_aux_updated = cone.dder3_aux_updated =
    cone.hess_fact_updated = false)

function setup_extra_data!(cone::EpiPerSepSpectral)
    @views cone.w_view = cone.point[3:end]
    setup_csqr_cache(cone)
    return cone
end

get_nu(cone::EpiPerSepSpectral) = 2 + cone.d

include("vectorcsqr.jl")
include("matrixcsqr.jl")

include("sepspectralfun.jl")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "xkbcommon"
version = v"0.9.1"

# Collection of sources required to build xkbcommon
sources = [
    ArchiveSource("https://xkbcommon.org/download/libxkbcommon-$(version).tar.xz",
                  "d4c6aabf0a5c1fc616f8a6a65c8a818c03773b9a87da9fbc434da5acd1199be0"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libxkbcommon-*/

# We need to run `wayland-scanner` on the host system
apk add wayland-dev

atomic_patch -p1 ../patches/meson_build.patch
atomic_patch -p1 ../patches/meson_options.patch

mkdir build && cd build
meson .. --cross-file="${MESON_TARGET_TOOLCHAIN}" \
    -Denable-docs=false \
    -Dnative-wayland-scanner="/usr/bin/wayland-scanner"
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [p for p in supported_platforms() if Sys.islinux(p)]

# The products that we will ensure are always built
products = [
    LibraryProduct("libxkbcommon", :libxkbcommon),
    LibraryProduct("libxkbcommon-x11", :libxkbcommon_x11),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"),
    Dependency("Xorg_xkeyboard_config_jll"),
    Dependency("Xorg_libxcb_jll"),
    Dependency("Wayland_jll"),
    Dependency("Wayland_protocols_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using TryCatch
using Test

@testset "TryCatch.jl" begin

    # Write your tests here.
    @test (@try 1 @else 2) == 2

    @test (@try @catch e @else 2) == 2

    @test (@try @catch @else 2) == 2

    @test (
        @try begin
            sqrt("0")
        @catch
            2
        end
    ) == 2

    @test (@try sqrt("0") @catch _ 2) == 2

    @test (@try 1) == 1

    @test_throws MethodError (@try sqrt("34"))

    @test_throws MethodError (@try sqrt("34") (@catch e::ErrorException e))

    # If we throw, the Exception block is returned
    @test (@try sqrt("34") @catch e::MethodError 5) == 5

    # catch MethodError
    @test (
        @try begin
            sqrt("34")
        @catch e isa MethodError
            52
        end
    ) == 52

    # catch MethodError
    @test (
        @try begin
            sqrt("34")
        @catch e::MethodError
            52
        end
    ) == 52

    # catch by complicated boolean expression
    @test (
        @try begin
            sqrt("34")
        @catch e->contains(string(e), "MethodError") 1
            52
            3
        end
    ) == 3

    # catch by complicated boolean expression
    @test (
    @try begin
        sqrt("34")
        @catch jasgdakja->contains(string(jasgdakja), "MethodError") begin
            52
        end
    end
    ) == 52

    # Let's see multiple catch blocks
    @test (
        @try begin
            sqrt("35")
        @catch _->false 15
        @catch _::ErrorException 20
        @catch _::MethodError 25
        @catch _->true 30
        end
    ) == 25

    # Finally annotation
    @test (begin
        a = 5
        @try begin
            @finally a = 6
        end
        a
    end) == 6

    # Finally annotation with catch
    @test (begin
        a = 5
        @try begin
            sqrt("35")

            @catch _::MethodError 25
            @finally a = 6
        end
        a
    end) == 6

    # Finally annotation with success annotation
    @test (begin

        a = nothing
        b = nothing

        @try begin
            sqrt("35")
        @catch _::MethodError
            25
        @else
            b = 9
        @finally
            a = 6
        end

        b
    end) === nothing

    # Finally annotation with success annotation
    @test (begin
        a = nothing
        b = nothing
        @try begin
        @catch _::MethodError
            25
        @else
            b = 9
        @finally
            a = 6
        end
        b
    end) === 9

    # Test syntax errors
    @test_throws LoadError @macroexpand @try 1 2 3
    @test_throws LoadError @macroexpand @try 3 @catch e 1 2
    @test @try true @catch e true @finally
    @test_throws LoadError @macroexpand @try true @finally @catch e true
    @test @try true @catch e true
    @test_throws LoadError @macroexpand @try true @catch 1+e true
    @test @try true @catch e::MethodError true

end


text:
# Calculate "G(tau, 0)", i.e. G(slice,1) as
#
# G(slice,1) = [B(slice, 1)^-1 + B(beta, slice)]^-1
#
# which is (in principle) equal to B(slice,1)G(1)
function calc_tdgf_udv(mc::AbstractDQMC, slice::Int, safe_mult::Int=mc.p.safe_mult)
  if slice != 1
    Ul, Dl, Vtl = calc_Bchain_inv_udv(mc, 1, slice-1, safe_mult)
  else
    Ul, Dl, Vtl = mc.s.eye_full, mc.s.ones_vec, mc.s.eye_full
  end

  if slice != mc.p.slices
    Ur, Dr, Tr = calc_Bchain_udv(mc, slice, mc.p.slices, safe_mult)
  else
    Ur, Dr, Tr = mc.s.eye_full, mc.s.ones_vec, mc.s.eye_full
  end

  # time displace
  U, D, Vt = inv_sum_udvs(Ul, Dl, Vtl, Ur, Dr, Tr)
  effective_greens2greens!(mc, U, Vt)

  rmul!(U, Diagonal(D))
  return U*Vt
end

global LEFT = true
global RIGHT = false
"""
Calculate UDVs at safe_mult time slices of
dir = LEFT:
inv=false:  B(tau, 1) = B(tau) * B(tau-1) * ... * B(1)                    # mult left, 1:tau
inv=true:   [B(tau, 1)]^-1 = B(1)^-1 * B(2)^-1 * ... B(tau)^-1            # mult inv right, 1:tau

udv[i] = from 1 to mc.s.ranges[i][end]

dir = RIGHT:
inv=false:  B(beta, tau) = B(beta) * B(beta-1) * ... * B(tau)             # mult right, beta:tau
inv=true:   [B(beta, tau)]^-1 = B(tau)^-1 * B(tau+1)^-1 * ... B(beta)^-1  # mult inv left, beta:tau

udv[i] = from mc.s.ranges[i][1] to mc.p.slices (beta)
"""
function calc_tdgf_B_udvs(mc::AbstractDQMC; inv::Bool=false, dir::Bool=LEFT)
  G = geltype(mc)
  flv = mc.p.flv
  N = mc.l.sites
  nranges= length(mc.s.ranges)
  curr_U_or_V = mc.s.curr_U
  eye_full = mc.s.eye_full
  ones_vec = mc.s.ones_vec
  ranges = mc.s.ranges

  u_stack = [zeros(G, flv*N, flv*N) for _ in 1:nranges]
  d_stack = [zeros(Float64, flv*N) for _ in 1:nranges]
  v_stack = [zeros(G, flv*N, flv*N) for _ in 1:nranges]

  rightmult = false
  ((dir == RIGHT && !inv) || (dir == LEFT && inv)) && (rightmult = true)

  # @show rightmult

  rng = 1:length(ranges)
  dir == RIGHT && (rng = reverse(rng))

  # Calculate udv[i], given udv[i-1]
  @inbounds for (i, rngidx) in enumerate(rng)

    if i == 1
      copyto!(curr_U_or_V, eye_full)
    else
      if !rightmult
        copyto!(curr_U_or_V, u_stack[i-1])
      else
        copyto!(curr_U_or_V, v_stack[i-1])
      end
    end

    slice_range = dir == RIGHT ? reverse(ranges[rngidx]) : ranges[rngidx]

    for slice in slice_range
      if inv == false
        if dir == LEFT
          multiply_B_left!(mc, slice, curr_U_or_V)
        else
          # rightmult
          multiply_B_right!(mc, slice, curr_U_or_V)
        end
      else
        if dir == LEFT
          # rightmult
          multiply_B_inv_right!(mc, slice, curr_U_or_V)
        else
          multiply_B_inv_left!(mc, slice, curr_U_or_V)
        end
      end
    end

    if i != 1
      if !rightmult
        rmul!(curr_U_or_V, Diagonal(d_stack[i-1]))
      else
        lmul!(Diagonal(d_stack[i-1]), curr_U_or_V)
      end
    end

    if !rightmult
      u_stack[i], d_stack[i], Vt = decompose_udv!(curr_U_or_V)
    else
      U, d_stack[i], v_stack[i] = decompose_udv!(curr_U_or_V)
    end

    if i == 1
      if !rightmult
        mul!(v_stack[i], Vt, eye_full)
      else
        mul!(u_stack[i], eye_full, U)
      end
    else
      if !rightmult
        mul!(v_stack[i], Vt, v_stack[i-1])
      else
        mul!(u_stack[i], u_stack[i-1], U)
      end
    end
  end

  if dir == LEFT
    return u_stack, d_stack, v_stack
  else
    return reverse(u_stack), reverse(d_stack), reverse(v_stack)
  end
end

function calc_tdgfs_udv!(mc)
  G = geltype(mc)
  M = mc.p.slices
  N = mc.l.sites
  flv = mc.p.flv
  Nflv = N * flv
  safe_mult = mc.p.safe_mult
  eye_full = mc.s.eye_full
  ones_vec = mc.s.ones_vec

  # allocate matrices if not yet done TODO: EVENTUALLY THIS SHOULD BE REMOVED
  try
    mc.s.meas.Gt0[1]
    mc.s.meas.G0t[1]
  catch
    allocate_tdgf!(mc)
  end

  Gt0 = mc.s.meas.Gt0
  G0t = mc.s.meas.G0t

  # ---- first, calculate Gt0 and G0t only at safe_mult slices
  # right mult (Gt0)
  BT0Inv_u_stack, BT0Inv_d_stack, BT0Inv_v_stack = calc_tdgf_B_udvs(mc, inv=true, dir=LEFT);
  BBetaT_u_stack, BBetaT_d_stack, BBetaT_v_stack = calc_tdgf_B_udvs(mc, inv=false, dir=RIGHT);

  # left mult (G0t)
  BT0_u_stack, BT0_d_stack, BT0_v_stack = calc_tdgf_B_udvs(mc, inv=false, dir=LEFT);
  BBetaTInv_u_stack, BBetaTInv_d_stack, BBetaTInv_v_stack = calc_tdgf_B_udvs(mc, inv=true, dir=RIGHT);

  safe_mult_taus = 1:safe_mult:mc.p.slices
  @inbounds for i in 1:length(safe_mult_taus) # i = ith safe mult time slice
    tau = safe_mult_taus[i] # tau = tauth (overall) time slice
    if i != 1
      # Gt0
      inv_sum_udvs!(mc, Gt0[tau], BT0Inv_u_stack[i-1], BT0Inv_d_stack[i-1], BT0Inv_v_stack[i-1],
                   BBetaT_u_stack[i], BBetaT_d_stack[i], BBetaT_v_stack[i]) # G(i,0) = G(mc.s.ranges[i][1], 0), i.e. G(21, 1) for i = 3
      effective_greens2greens!(mc, Gt0[tau])

      # G0t
      inv_sum_udvs!(mc, G0t[tau], BT0_u_stack[i-1], BT0_d_stack[i-1], BT0_v_stack[i-1],
                   BBetaTInv_u_stack[i], BBetaTInv_d_stack[i], BBetaTInv_v_stack[i]) # G(i,0) = G(mc.s.ranges[i][1], 0), i.e. G(21, 1) for i = 3
      effective_greens2greens!(mc, G0t[tau])
    else
      # Gt0
      inv_one_plus_udv_scalettar!(mc, Gt0[tau], BBetaT_u_stack[1], BBetaT_d_stack[1], BBetaT_v_stack[1])
      effective_greens2greens!(mc, Gt0[tau])

      # G0t
      inv_one_plus_udv_scalettar!(mc, G0t[tau], BBetaTInv_u_stack[1], BBetaTInv_d_stack[1], BBetaTInv_v_stack[1])
      effective_greens2greens!(mc, G0t[tau]) # TODO: check analytically that we can still do this
    end
  end

  # ---- fill time slices between safe_mult slices
  fill_tdgf!(mc, Gt0, G0t)

  @inbounds for i in 1:M
    G0t[i] .*= -1 # minus sign
  end

  nothing
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
const _hermite = OrthogonalPolynomialsQuasi.Hermite()
# perhaps call fastGaussQuadrature directly
const DEFAULT_HERMITE_INTEGRATOR = sqrt(2π)*expectation(Normal(); n=100)
const DEFAULT_SQUARED_HERMITE_INTEGRATOR = sqrt(π)*expectation(Normal(0, 1/sqrt(2)); n=101)

# physicist's hermite polynomials

function hermite_fun_noexp(xs, n)
	xs = big.(xs)
	n = big.(n )
	Float64.(_hermite[xs, n+1]./sqrt(2^n*factorial(n)*sqrt(pi)))
end

function hermite_fun(xs, n)
	hermite_fun_noexp(xs, n).*exp.(-xs.^2/2)
end

hermite_cf_coefficient(n) = im^n*sqrt(2*π)

#fourier_transform(t, order) = hermite_fun(t, order)*hermite_cf_coefficient(order)

# Hermite = Mixture of Hermite Basis Functions
struct HermiteBasisFunction
	q::Integer
end

(herm::HermiteBasisFunction)(x) = MCEB.hermite_fun(x, herm.q)

function expectation(herm::HermiteBasisFunction,
	                 base_integrator::Expectation = DEFAULT_HERMITE_INTEGRATOR)
	hermite_integrator = deepcopy(base_integrator)
	wts = hermite_fun_noexp.(nodes(hermite_integrator), herm.q)
	Expectations.weights(hermite_integrator) .*= wts
	hermite_integrator
end

# wow this should really be dispatching on integrator
function (target::MarginalDensityTarget{<:DiscretizedStandardNormalSample})(integrator::Expectation)
    Z_disc = location(target)
    grid = Z_disc.samples.mhist.grid
    i = Z_disc.bin
    if i == 1
        base_integrand = μ -> cdf(Normal(μ), first(grid))
    elseif i == lastindex(Z_disc.samples.mhist)
        base_integrand = μ -> ccdf(Normal(μ), last(grid))
    else
        base_integrand = μ -> (cdf(Normal(μ), grid[i]) -  cdf(Normal(μ), grid[i-1]))
    end
    # refactor into riesz representer.. -> Oh already thought about it...
    integrator(base_integrand)
end

# Here ``g^*(t) = \\int \\exp(2 i \\pi \\mu) g(\\mu)d\\mu$ is the Fourier transform of ``g``.
"""
	HermitePriorClass(qmax, sobolev_order, sobolev_bound, solver)

Class of densities ``g \\in \\mathcal{G}`` that can be written as ``g(\\mu) = \\sum_{j=0}^{qmax} h_j(\\mu)``, where
``h_j`` is the ``j``-th Hermite function and further satisfy the Sobolev constraint
(with `sobolev_order` corresponding to ``b`` and `sobolev_bound` to ``C``):

```math
\\int_{-\\infty}^{\\infty} |g^*(t)|^2(t^2+1)^b dt \\leq 2\\pi C
```
The `solver` object is a `MathOptInterface` compatible optimizer such as `Mosek.Optimizer` that will be
used to solve the convex modulus of continuity (or worst-case-bias estimation) problem.
"""
Base.@kwdef mutable struct HermitePriorClass{T<:Real} <: ConvexPriorClass
	qmax::Integer
	sobolev_order::T  #not necessarily but should suffice for now
	sobolev_bound::Float64
	integrators = expectation.(HermiteBasisFunction.(0:qmax), Ref(DEFAULT_HERMITE_INTEGRATOR))
	squared_integrator = DEFAULT_SQUARED_HERMITE_INTEGRATOR
	solver = nothing
	solver_params = []
end

function add_prior_variables!(model::JuMP.Model, hermite_class::HermitePriorClass; var_name = "πs") # adds constraints
    qmax = hermite_class.qmax
	nparams = qmax + 1
    tmp_vars = @variable(model, [i=1:nparams])
    model[Symbol(var_name)] = tmp_vars
	cf_coefs = hermite_cf_coefficient.(0:qmax) .* hermite_fun_noexp.(0.0, 0:qmax)
	con = @constraint(model, dot(tmp_vars, real.(cf_coefs)) == 1.0)
	con2 = @constraint(model, dot(tmp_vars, imag.(cf_coefs)) == 0.0)
	#tmp_grid = -7:0.02:4
	tmp_grid = sort([nodes(hermite_class.squared_integrator); nodes(hermite_class.integrators[1])])
	# rewrite below
	A_mat_prior = [hermite_fun(x, qidx)  for x in tmp_grid, qidx=0:qmax]
	nonneg_constraint = @constraint(model, A_mat_prior*tmp_vars .>= 0)

	integrator2 = hermite_class.squared_integrator
	b = hermite_class.sobolev_order
	sobolev_bound = sqrt(hermite_class.sobolev_bound*2π)

	ft_mat = [hermite_fun_noexp(x, qidx)*hermite_cf_coefficient(qidx)  for x in integrator2.nodes, qidx=0:qmax]
	real_ft = (real(ft_mat) * tmp_vars) .* sqrt.(integrator2.weights .* (integrator2.nodes.^2 .+ 1).^b)
	imag_ft = imag(ft_mat) * tmp_vars .* sqrt.(integrator2.weights .* (integrator2.nodes.^2 .+ 1).^b)

	@constraint(model, [sobolev_bound; real_ft; imag_ft] in SecondOrderCone())

	#@constraint(model, sum( real_ft.^2 .* integrator2.weights) + sum( imag_ft.^2 .* integrator2.weights) <= sobolev_bound)

    tmp_vars
end

function marginalize(hermite_class::HermitePriorClass,
                     Zs_discr::DiscretizedStandardNormalSamples,
                     param_vec) #-> A*param_vec
	qdtrs =  hermite_class.integrators
	Zs_discr_length = length(Zs_discr.mhist.grid) + 1
	A_mat = [ MarginalDensityTarget(Zs_discr(i))(qdtr) for i=1:Zs_discr_length, qdtr in qdtrs]
    A_mat*param_vec
end

function linear_functional(hermite_class::HermitePriorClass,
                           target::PriorDensityTarget, #let's hardcode to this for now.
                           param_vec)
	x = location(target)
	qmax = hermite_class.qmax
    v = hermite_fun.(x, 0:qmax)
    dot(v, param_vec)
end

struct HermiteQuasiDistribution<: Distribution{Univariate, Continuous}
	q::Integer
	coefs::Vector{Float64}
end

function pdf(herm_dbn::HermiteQuasiDistribution, x)
	 dot(hermite_fun.(x, 0:herm_dbn.q), herm_dbn.coefs)
end

function cf(herm_dbn::HermiteQuasiDistribution, t)
	 cf_coefs = hermite_cf_coefficient.(0:herm_dbn.q)
	 dot(hermite_fun.(t, 0:herm_dbn.q), herm_dbn.coefs .* cf_coefs)
end

function marginalize(hermite_dbn::HermiteQuasiDistribution,
                     Zs_discr::DiscretizedStandardNormalSamples)
	#temporary hack
   hm_prior_class = HermitePriorClass(qmax=hermite_dbn.q,
                    sobolev_order =2 , sobolev_bound = 10.0)
   ws = marginalize(hm_prior_class, Zs_discr, hermite_dbn.coefs)
   mhist = Zs_discr.mhist
   mhist = @set mhist.hist.weights = ws
end

function (hermite_class::HermitePriorClass)(param_vec)
    q = hermite_class.qmax
  	HermiteQuasiDistribution(q, param_vec)
end

# temporary hacks
struct GaussianHermite
	f
end

function marginalize(hermite_dbn::HermiteQuasiDistribution, Z::StandardNormalSample)
  αs = hermite_dbn.coefs
  qmax = hermite_dbn.q
  GaussianHermite(x -> sum([αs[i+1]*_tmp_conv(i,x) for i in 0:qmax]))
end

function _tmp_conv(j, x)
	j = big(j)
	big(1/sqrt(2*sqrt(pi)*2^j*factorial(j))*exp(-x^2/4)*x^j)
end
pdf(g::GaussianHermite, x) = Float64(g.f(x))
broadcastable(g::GaussianHermite) = Ref(g)


text:
#load pendulum and controllers
include("../FiniteDifferences/heavy_chain.jl")
include("feedforward.jl")
include("feedback.jl")

#define parameters
enable_anim = false
enable_sim = true
#x_coords = collect(linspace(0, 1.178, 200))
x_coords = collect(1.0-(logspace(1, 0, 200)-1)./9)*1.178
chain = Chain(x_coords)
dt = 1/2000
all_t = collect(0:dt:15)
(F_ff, xL_target, x0_target, v0_target, phi0_target, phi0_dot_target) = chain_feedforward(chain, all_t)

#plot signals
using Plots
pyplot()
fig = plot(layout = (3, 1))
plot!(fig[1], all_t, xL_target, lab = "xL")
plot!(fig[1], all_t, x0_target, l = :dot, lab = "x0")
plot!(fig[2], all_t, phi0_target, lab = "phi")
plot!(fig[2], all_t, phi0_dot_target, l = :dot, lab = "dphidt")
plot!(fig[3], all_t, F_ff, lab = "F")
gui()

#simulate
if enable_sim
  println("Starting simulation")
  all_xL = Array{Float64}(length(xL_target))
  all_x0 = Array{Float64}(length(x0_target))
  for i = 1:length(all_t)
    u_in = F_ff[i]
    all_x0[i] = chain.w[1]
    all_xL[i] = chain.w[end]
    simulate!(chain, u_in, dt)
  end
  fig = plot(layout = (2, 1))
  plot!(fig[1], all_t, x0_target, l=:dot, lab = "x0_target")
  plot!(fig[1], all_t, all_x0, lab = "x0", c=:black)

  plot!(fig[2], all_t, xL_target, l=:dot, lab = "xL_target")
  plot!(fig[2], all_t, all_xL, lab = "xL", c=:black)
end

#animate
if enable_anim
  println("Starting animation")
  include("../heavy_chain_plot.jl")
  xmin = -0.2
  xmax = 0.7
  chain = Chain(x_coords)
  anim = Animation()
  all_u = [NaN for i = 1:length(all_t)]
  iter=1
  u_in = 1.0
  for t = all_t
    #take frame
    if mod(t, 0.04)<dt*0.5 || mod(t, 0.04)>(0.04-dt*0.5)
      println("Frame at t = $(t)")
      l = Plots.@layout [a{0.7h}; b]
      fig = plot(layout = l)
      heavy_chain_plot!(fig[1], chain, xmin, xmax, t)
      plot!(fig[2], all_t, all_u, lab = "u")
      frame(anim)
    end

    #simulate one step
    u_in = F_ff[iter]
    simulate!(chain, u_in, dt)
    all_u[iter] = u_in
    iter+=1
  end
  gif(anim, "animation.mp4", fps = round(Int, 1/0.04))
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Ray of light propagating in space for the raytracing algorithms

"""
    Ray

A ray of light propagating in space.

# Fields

- `origin::Point`: the ([`Point`](@ref)) where the ray originated.
- `dir::Vec`: a ([`Vec`](@ref)) representing the direction along which this ray propagates.
- `tmin::Float32`: the minimum distance travelled by the ray is this number times `dir`.
- `tmax::Float32`: the maximum distance travelled by the ray is this number times `dir`.
- `depth::Int`: number of times this ray was reflected/refracted.

See also: [`Ray(::Float32)`](@ref)
"""
struct Ray
    origin::Point
    dir::Vec
    tmin::Float32
    tmax::Float32
    depth::Int

    function Ray(origin::Point, dir::Vec, tmin::Float32, tmax::Float32, depth::Int)
        tmin < tmax || throw(ArgumentError("`tmin >= tmax` is not allowed: `tmin = $tmin`, `tmax = $tmax`"))
        new(origin, dir, tmin, tmax, depth)
    end
end

@doc """
    Ray(origin::Point, dir::Vec, tmin::Float32, tmax::Float32, depth::Int)

Constructor for a [`Ray`](@ref) instance.
""" Ray(::Point, ::Vec, ::Float32, ::Float32, ::Int)

"""
    Ray(origin::Point, dir::Vec
        ; tmin::Float32 = 1f-5,
          tmax::Float32 = Inf32,
          depth::Int = 0)

Constructor for a [`Ray`](@ref) instance.
"""
Ray(origin::Point, dir::Vec; tmin::Float32 = 1f-5, tmax::Float32 = Inf32, depth::Int = 0) = Ray(origin, dir, tmin, tmax, depth)

@doc raw"""
    (r::Ray)(t::Float32)

Return a [`Point`](@ref) lying on the given [`Ray`](@ref) at `t`.

An instance of `Ray` can be called as a function returning a `Point` given the position parameter `t`:

```math
\mathrm{ray\_origin} + \mathrm{ray\_direction} \cdot t
```

Argument `t` must be included between `r.tmin` and `r.tmax` or be equal to 0.
If `t` is zero, then the returned point is the origin of `r`.

# Examples

```jldoctest
julia> ray = Ray(ORIGIN, VEC_X)
Ray
 ↳ origin = Point(0.0, 0.0, 0.0)
 ↳ dir    = Vec(1.0, 0.0, 0.0)
 ↳ tmin   = 1.0e-5
 ↳ tmax   = Inf
 ↳ depth  = 0

julia> ray(5f0)
Point with eltype Float32
x = 5.0, y = 0.0, z = 0.0
```
"""
function (r::Ray)(t::Float32)
    t == 0 && return r.origin
    r.tmin <= t <= r.tmax || throw(ArgumentError("argument `t` must have a value between `r.tmin = $(r.tmin)` and `r.tmax = $(r.tmax)`: got `t = $t`"))
    r.origin + r.dir * t
end

################
# Miscellaneous

function show(io::IO, ::MIME"text/plain", ray::T) where {T <: Ray}
    print(io, T)
    fns = fieldnames(T)
    n = maximum(fns .|> String .|> length)
    for fieldname ∈ fns
        println(io)
        print(io, " ↳ ", rpad(fieldname, n), " = ", getfield(ray, fieldname))
    end
end

eltype(::Ray) = Float32
eltype(::Type{Ray}) = Float32

"""
    ≈(r1::Ray, r2::Ray)

Check if two [`Ray`](@ref) represent the same ray of light or not.
"""
(≈)(r1::Ray, r2::Ray) = (r1.origin ≈ r2.origin) && (r1.dir ≈ r2.dir)

"""
    *(t::Transformation, r::Ray)

Transform a [`Ray`](@ref) with the given [`Transformation`](@ref).

# Examples

```jldoctest
julia> ray = Ray(ORIGIN, VEC_X)
Ray
 ↳ origin = Point(0.0, 0.0, 0.0)
 ↳ dir    = Vec(1.0, 0.0, 0.0)
 ↳ tmin   = 1.0e-5
 ↳ tmax   = Inf
 ↳ depth  = 0

julia> Transformation() * ray
Ray
 ↳ origin = Point(0.0, 0.0, 0.0)
 ↳ dir    = Vec(1.0, 0.0, 0.0)
 ↳ tmin   = 1.0e-5
 ↳ tmax   = Inf
 ↳ depth  = 0

julia> translation(2,4,-6) * ray
Ray
 ↳ origin = Point(2.0, 4.0, -6.0)
 ↳ dir    = Vec(1.0, 0.0, 0.0)
 ↳ tmin   = 1.0e-5
 ↳ tmax   = Inf
 ↳ depth  = 0
```
"""
(*)(t::Transformation, r::Ray) = Ray(t*r.origin, t*r.dir, r.tmin, r.tmax, r.depth)


text:
baremodule Atomix

macro atomic end
macro atomicreplace end
macro atomicswap end

abstract type IndexableRef{Indexable} end

function asstorable end

# Maybe it's useful to make `pointer` customizable for those who want to
# dispatch at the level of UnsafeAtomics?
function pointer end

function gcroot end

function get end
function set! end
function modify! end
function swap! end
function replace! end

module Internal

import ..Atomix: @atomic, @atomicswap, @atomicreplace
using ..Atomix

using Base.Meta: isexpr
using Base: @propagate_inbounds
using UnsafeAtomics:
    Ordering, UnsafeAtomics, monotonic, acquire, release, acq_rel, seq_cst, right

include("utils.jl")
include("references.jl")
include("generic.jl")
include("core.jl")
include("sugar.jl")

function define_docstring()
    path = joinpath(@__DIR__, "..", "README.md")
    include_dependency(path)
    doc = read(path, String)
    doc = replace(doc, r"^```julia"m => "```jldoctest README")
    # Setting `LineNumberNode` to workaround an error from logging(?) `no method
    # matching getindex(::Nothing, ::Int64)`:
    ex = :($Base.@doc $doc Atomix)
    ex.args[2]::LineNumberNode
    ex.args[2] = LineNumberNode(1, Symbol(path))
    Base.eval(Atomix, ex)
end

end  # module Internal

const Ordering = Internal.Ordering

const monotonic = Internal.monotonic
const acquire = Internal.acquire
const release = Internal.release
const acq_rel = Internal.acq_rel
const seq_cst = Internal.seq_cst

# Julia names
const acquire_release = acq_rel
const sequentially_consistent = seq_cst

const right = Internal.UnsafeAtomics.right

Internal.define_docstring()

end  # baremodule Atomix


text:
module MassSpec

import CSV
import Statistics

using DataFrames
using Printf
using RecipesBase
using StatsBase

struct MassSpectrum{T<:Number} <: AbstractArray{Tuple{T,T},1}
    mass::Array{T,1}
    intensity::Array{T,1}
end

intensity(ms::MassSpectrum) = maximum(ms.intensity)
Base.sum(ms::MassSpectrum) = sum(ms.intensity)
Base.size(ms::MassSpectrum) = (size(ms.mass)[1],)
Base.getindex(ms::MassSpectrum, i::Union{Integer, CartesianIndex}) = (ms.mass[i], ms.intensity[i])
Base.getindex(ms::MassSpectrum, i) = MassSpectrum(ms.mass[i], ms.intensity[i])

function Base.setindex!(ms::MassSpectrum, val, i::Integer)
    setindex!(ms.mass, val[1], i)
    setindex!(ms.intensity, val[2], i)
end

function Base.setindex!(ms::MassSpectrum, val, i::Union{AbstractRange, Colon})
    setindex!(ms.mass, val[:,1], i)
    setindex!(ms.intensity, val[:,2], i)
end

function Base.deleteat!(ms::MassSpectrum, x)
    deleteat!(ms.mass, x)
    deleteat!(ms.intensity, x)
end

function MassSpectrum(path::AbstractString)
    df = CSV.read(path, header=["mass", "mV"], datarow=2, delim=',', allowmissing=:none)
    MassSpectrum(df.mass, df.mV)
end

mass_format(m) = Printf.@sprintf "%.1f" m

include("operations.jl")
include("plotting.jl")

export MassSpectrum, Mass, discretize, undiscretize, subtract, distance

end


text:
import CLBLAS
import OpenCL
    const clblas = CLBLAS
    clblas.setup()
    using OpenCL: cl

    device, ctx, queue = clblas.get_next_compute_context()

    N = Unsigned(5)
    alpha1 = cl.cl_float(5)
    alpha2 = cl.cl_float(2)
    data = Vector{Float32}(5)
    data[1] = 0.1
    data[2] = 0.22
    data[3] = 0.333
    data[4] = 0.4444
    data[5] = 0.55555

    X = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=data)
    offx = Unsigned(0)
    incx = cl.cl_int(1)
    ncq = cl.cl_uint(1)
    clq = queue.id
    newl = cl.cl_uint(0)

    ptrEvent = clblas.clblasSscal(N, alpha1, X.id, offx, incx, [queue])

    ptrEvent2 = clblas.clblasSscal(N, alpha2, X.id, offx, incx,[queue], ptrEvent)

    cl.api.clWaitForEvents(cl.cl_uint(1), ptrEvent2)

    result = Vector{Float32}(length(X))
    cl.enqueue_read_buffer(queue, X, result, Unsigned(0), nothing, true)

    expected = data * alpha1 * alpha2
    println(result)
    println("------------")
    println(expected)

    try
        @assert(expected == result)
        info("success!")
    finally
        clblas.teardown()
    end


text:
using Test
using DataWrangler: normality

@testset "normality" begin

    x = 1:100
    @test normality(x) ≈ 6.002400480071892

    x = randn(100)
    @test normality(x) < 6

end


text:
"""Class for holding IR bases and sparse-sampling objects

An object of this class holds IR bases for fermion and bosons
and associated sparse-sampling objects

Attributes:
    basis_f (FiniteTempBasis):
        Fermion basis
    basis_b (FiniteTempBasis):
        Boson basis
    beta (Float64):
        Inverse temperature
    wmax (Float64):
        Cut-off frequency
    tau (Vector{Float64}):
        Sampling points in the imaginary-time domain
    wn_f (Vector{Int64}):
        Sampling fermionic frequencies
    wn_b (Vector{Int64}):
        Sampling bosonic frequencies
    smpl_tau_f (TauSampling):
        Sparse sampling for tau & fermion
    smpl_tau_b (TauSampling):
        Sparse sampling for tau & boson
    smpl_wn_f (MatsubaraSampling):
        Sparse sampling for Matsubara frequency & fermion
    smpl_wn_b (MatsubaraSampling):
        Sparse sampling for Matsubara frequency & boson
    sve_result (Tuple{PiecewiseLegendrePoly,Vector{Float64},PiecewiseLegendrePoly}:
        Results of SVE
"""
struct FiniteTempBasisSet
    basis_f::FiniteTempBasis
    basis_b::FiniteTempBasis
    smpl_tau_f::TauSampling
    smpl_tau_b::TauSampling
    smpl_wn_f::MatsubaraSampling
    smpl_wn_b::MatsubaraSampling
end

"""
Create basis sets for fermion and boson and
associated sampling objects.
Fermion and bosonic bases are constructed by SVE of the logistic kernel.
"""
function FiniteTempBasisSet(beta::Real, wmax::Real, eps::Real;
                            sve_result::Union{SVEResultType,Nothing}=nothing)
    if sve_result === nothing
        # Create bases by sve of the logistic kernel
        basis_f, basis_b = finite_temp_bases(beta, wmax, eps)
    else
        # Create bases using the given sve results
        basis_f = FiniteTempBasis(fermion, beta, wmax, eps; sve_result=sve_result)
        basis_b = FiniteTempBasis(boson, beta, wmax, eps; sve_result=sve_result)
    end

    return FiniteTempBasisSet(basis_f, basis_b,
                              TauSampling(basis_f), TauSampling(basis_b),
                              MatsubaraSampling(basis_f), MatsubaraSampling(basis_b))
end

function Base.getproperty(bset::FiniteTempBasisSet, d::Symbol)
    if d === :beta
        return getfield(bset, :basis_f).beta
    elseif d === :wmax
        return getfield(bset, :basis_f).wmax
    elseif d === :tau
        return getfield(bset, :smpl_tau_f).sampling_points
    elseif d === :wn_f
        return getfield(bset, :smpl_wn_f).sampling_points
    elseif d === :wn_b
        return getfield(bset, :smpl_wn_b).sampling_points
    elseif d === :sve_result
        return getfield(bset, :basis_f).sve_result
    else
        return getfield(bset, d)
    end
end

function Base.propertynames(bset::FiniteTempBasisSet, private::Bool=false)
    return (:beta,
            :wmax,
            :tau,
            :wn_f,
            :wn_b,
            :sve_result,
            fieldnames(FiniteTempBasisSet)...)
end


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

# # Assembling stiffness matrix and force vector for 2d heat problem

using FEMBase, HeatTransfer, Test, SparseArrays

X = Dict(
    1 => [0.0,0.0],
    2 => [1.0,0.0],
    3 => [1.0,1.0],
    4 => [0.0,1.0])

T = Dict(
    1 => 0.0,
    2 => 0.0,
    3 => 0.0,
    4 => 0.0)

element = Element(Quad4, [1, 2, 3, 4])
update!(element, "geometry", X)
update!(element, "temperature", 0.0 => T)
update!(element, "thermal conductivity", 6.0)
update!(element, "heat source", 4.0)

problem = Problem(PlaneHeat, "test problem", 1)
add_elements!(problem, [element])
assemble!(problem, 0.0)
K = Matrix(sparse(problem.assembly.K))
f = Vector(sparse(problem.assembly.f)[:])
K_expected = [
               4.0 -1.0 -2.0 -1.0
              -1.0  4.0 -1.0 -2.0
              -2.0 -1.0  4.0 -1.0
              -1.0 -2.0 -1.0  4.0
             ]
f_expected = [1.0; 1.0; 1.0; 1.0]
@test isapprox(K, K_expected)
@test isapprox(f, f_expected)
@test get_unknown_field_name(problem) == "temperature"


text:
function thin(proj_massfunc, θ, rs, tgα, N=100; ext=1.0)
    const cθ = cos(θ)
    const sθ = sin(θ)
    const x0 = -rs*cθ
    const y0 = rs*sθ
    const b = y0-tgα*x0
    const α=atan(tgα)
    const M0p = proj_massfunc(b)
    const deflection_angle = 4*M0p*G/C/C/b
    const nα = α-deflection_angle
    const tgnα = tan(nα)
    nx = -b/tgnα
    if nα >0.0
        nx = rs
    end

    lx = linspace(0.0,1.0,div(N,2)).^3
    x = [reverse(lx).*x0; lx.*nx.*ext]
    y = zeros(length(x))
    for i in 1:length(x)
        y[i] = if x[i]<0.0
            y0+tgα*(x[i]-x0)
        else
            b+tgnα*x[i]
        end
    end

    in_angle =  atan(tgα)
    out_angle =  nα

    u = 1./(x.*x+y.*x)
    θ = atan2(y, x)
    du = dudθ(θ, tgα, 1./u) # TODO: test this du is probably wrong, and not currently used

    x, y, u, du, in_angle, out_angle, deflection_angle, π-atan2(y,x)
end

# some convenience functions:

thin_tiso200b(z, M200b, θ, rs, tgα, N=100; ext=1.0) =
    thin(r->tiso_mp(M200b, r, R200b(z, M200b)), θ, rs, tgα, N; ext=ext)
thin_tiso200c(z, M200c, θ, rs, tgα, N=100; ext=1.0) =
    thin(r->tiso_mp(M200c, r, R200c(z, M200c)), θ, rs, tgα, N; ext=ext)
thin_tisoVir(z, Mvir, θ, rs, tgα, N=100; ext=1.0) =
    thin(r->tiso_mp(Mvir, r, RVIR(z, MVir)), θ, rs, tgα, N; ext=ext)

thin_tnfw200b(z, M200b, Rscale, θ, rs, tgα, N=100; ext=1.0) =
    thin(r->tnfw_mp(M200b, Rscale, r, R200b(z, M200b)), θ, rs, tgα, N; ext=ext)
thin_tnfw200c(z, M200c, Rscale, θ, rs, tgα, N=100; ext=1.0) =
    thin(r->tnfw_mp(M200c, Rscale, r, R200c(z, M200c)), θ, rs, tgα, N; ext=ext)
thin_tnfwVir(z, Mvir, Rscale, θ, rs, tgα, N=100; ext=1.0) =
    thin(r->tnfw_mp(Mvir, Rscale, r, RVIR(z, Mvir)), θ, rs, tgα, N; ext=ext)

max_stars_repo_name: skariel/Lensing

max_issues_repo_name: skariel/Lensing

max_forks_repo_name: skariel/Lensing


text:
"""
    plot_missing(data::PopObj)
Return an interactive plot of the number of missing loci in individuals of a
`PopObj`, along with the number of missing individuals per locus.
"""
function plot_missing(data::PopObj)
    by_sample,by_loci = missing(data);
    #ys = Array[subdf[!, :missing] for subdf in groupby(by_sample[!, 1:3], :population)]
    #texts = Array[subdf[!, :name] for subdf in groupby(by_sample[!, 1:3], :population)]
    popnum = length(by_sample.population |> unique)
    colors = ["hsl($i, 50%, 50%)" for i in range(0, stop=300, length=popnum)]

    byind = box(
        by_sample,
        x = :missing,
        group = :population,
        marker_color=colors,
        text = :name,
        jitter = 1,
        pointpos = 0,
        marker_size = 8.5,
        boxpoints = "all",
        marker=attr(line=attr(width=0.75))
    )

    layout_ind = Layout(
        title = "Number of missing loci per population",
        hovermode = "closest",
        yaxis = attr(title = "# missing loci", zeroline = false),
        xaxis = attr(title = "Population", zeroline = false)
    )

    loci_hist = histogram(x = by_loci[!, :missing],
                          marker_color = "rgb(217, 217, 217)",
                          name = "",
                          text = "loci",
                          xbins_size = 1,
                          showlegend = false)

    layout_hist = Layout(title = "Distribution of Locus Missingness",
                         xaxis = attr(title = "# of times locus missing", zeroline = false),
                         yaxis = attr(title = "# of loci", zeroline = false)
                         )
    ind_plot =  plot(byind, layout_ind)
    loci_plot = plot(loci_hist, layout_hist)

    return [ind_plot loci_plot]
end

"""
    plot_locations(x::PopObj; region::String = "world", projection::String = "orthographic")

Returns a simple lower resolution interactive scatterplot of the individuals in a
`PopObj`. Default `region` and `projection` are "world" and "orthographic",
respectively. If a specific region is set, the plot will default to "mercator"
projection unless `projection =` is used to specify a different one.

Example:\n
`plot_locations(manatees, region = "usa", projection = "albers usa")`

[regions]\n
"usa","europe", "asia", "africa", "north america", "south america"

[projections]\n
"equirectangular", "mercator", "orthographic", "natural earth",
"kavrayskiy7", "miller", "robinson", "eckert4",
"azimuthal equal area", "azimuthal equidistant",
"conic equal area", "conic conformal" , "conic equidistant",
"gnomonic", "stereographic", "mollweide", "hammer",
"transverse mercator", "albers usa", "winkel tripel",
"aitoff", "sinusoidal"
"""
function plot_locations(x::PopObj; region::String = "world", projection::String = "orthographic")
    # test for missing?
    if projection ∉ ["equirectangular", "mercator", "orthographic", "natural earth",
                    "kavrayskiy7", "miller", "robinson", "eckert4",
                    "azimuthal equal area", "azimuthal equidistant",
                    "conic equal area", "conic conformal" , "conic equidistant",
                    "gnomonic", "stereographic", "mollweide", "hammer",
                    "transverse mercator", "albers usa", "winkel tripel",
                    "aitoff", "sinusoidal"]
        error("Projection not recognized. Please see the help doc for list of projection options")
    end
    y = PopOpt(x)
    popnum = length(y.samples.population |> unique)
    colors = ["hsl($i, 50%, 50%)" for i in range(0, stop=300, length=popnum)]
    df_split = groupby(y.samples, :population)
    map_scatter = [scattergeo(lat=df_split[i][!, :latitude],
                       lon=df_split[i][!, :longitude],
                       marker_line_color="rgb(62,90,112)", marker_line_width=1,
                       marker_color = colors[i],
                       name = df_split[i][!, :population][1],
                       text = df_split[i][!, :name]
                       ) for i in 1:popnum]
   if region == "world"
       geo = attr(scope = region,
                  projection_type = projection,
                  showcoastlines = false,
                  showcountries = true,
                  countrywidth = 0.75,
                  countrycolor = "rgb(255,255,255)",
                  subunitcolor = "rgb(255,255,255)",
                  showland = true,
                  landcolor = "rgb(217, 217, 217)",
                  )
   elseif lowercase(region) ∈ ["usa","europe", "asia", "africa", "north america", "south america"]
      if projection == "orthographic"
          proj = "mercator"
      else
          proj = projection
      end
      geo = attr(scope = lowercase(region),
                 showlakes = true,
                 lakecolor = "#a5a5a5",
                 showrivers = true,
                 rivercolor = "#a5a5a5",
                 showsubunits = true,
                 projection_type = proj,
                 showcoastlines = false,
                 showcountries = true,
                 countrywidth = 0.75,
                 countrycolor = "rgb(255,255,255)",
                 subunitcolor = "rgb(255,255,255)",
                 showland = true,
                 landcolor = "rgb(217, 217, 217)",
                 )
   else
       error("Leave \"region =\" empty or use one of: usa, europe, asia, africa, north america, south america")
   end
    layout = Layout(;title="Sample locations", showlegend=true, geo=geo)
    plot(map_scatter, layout)
end


text:
function gensym_index(column_names::NTuple{N, ColumnName}) where N
    safe_column_names = ntuple(i -> gensym(), length(column_names))
    mapping = Dict{ColumnName, Symbol}(column_names .=> safe_column_names)
    safe_column_names, mapping
end


text:
for standardize in (true, false)
    for bias in (0, 1)
        cntr = center(data; standardize=standardize, bias=bias)[1]
        wcntr = center(data, uniform_weights; standardize=standardize, bias=bias)[1]
        @test tolerance(cntr, wcntr)
    end
end


text:
"""
Adapted from sandreza/Learning/sandbox/gaussian_process.jl
https://github.com/sandreza/Learning/blob/master/sandbox/gaussian_process.jl
Changed handling of kernel functions; changed some variable names;
added log marginal likelihood function.
"""

using LinearAlgebra

"""
GP_multiple
# Description
- data structure for GPR computations where each gridpoint in the prediction has a different predictor
# Data Structure and Description
    GPs, Array of GP objects
    kernel, Kernel object
    x_train
"""
struct GP_multiple
    GPs::Array{GP}
    kernel::Kernel
    x_train
end

function stencil_range(D, stencil_size, i)
    k = Int(floor(stencil_size/2))
    if i-k < 1
        return 1:stencil_size
    elseif i-k+stencil_size-1 > D
        return D-stencil_size+1:D
    else
        start = i-k
        start:start+stencil_size-1
    end
end

stencil(stencil_range, data) = [x[stencil_range] for x in data]

function model(𝒟::ProfileData; kernel::Kernel = Kernel(), stencil_size=)
    # create instance of GP using data from ProfileData object
    stencil_ranges = [stencil_range(D,stencil_size,i) for i=1:D]
    GPs = [model(stencil(𝒟.x_train,r), stencil(𝒟.y_train,r), kernel, 𝒟.zavg[r]) for range in stencil_ranges]
    return GP_multiple(GPs, kernel, 𝒟.x_train);
end

"""
model_output(x, 𝒢::GP_multiple)
# Description
- Given state x, GP 𝒢, returns the mean GP prediction
# Arguments
- `x`: single scaled state
- `𝒢`: GP_multiple object with which to make the prediction
# Return
- `y`: scaled prediction
"""
function model_output(x, 𝒢::GP_multiple)
    return [model_output(x[i],𝒢.GPs[i]) for i in 1:length(x)]
end

"""
uncertainty(x, 𝒢::GP)
# Description
- Given state x and GP 𝒢, output the variance at a point
# Arguments
- `x`: state
# Return
- `var`: variance
"""
function uncertainty(x, 𝒢::GP)
    tmpv = zeros(size(𝒢.x_train)[1])
    for i in eachindex(𝒢.x_train)
        tmpv[i] = 𝒢.kernel(x, 𝒢.x_train[i])
    end
    # no ldiv for suitesparse
    tmpv2 = 𝒢.CK \ tmpv
    var = k(x, x) .- tmpv'*tmpv2  # var(f*) = k(x*,x*) - tmpv'*tmpv2
    return var
end

"""
compute_kernel_matrix(kernel, x)
# Description
- Computes the kernel matrix for GPR
# Arguments
- `k` : (Kernel) kernel function k(a,b).
- `x` : (array of predictors). x[1] is a vector
# Return
- `sK`: (symmetric matrix). A symmetric matrix with entries sK[i,j] = k(x[i], x[j]). This is only meaningful if k(x,y) = k(y,x) (it should)
"""
function compute_kernel_matrix(k, x)

    K = [k(x[i], x[j]) for i in eachindex(x), j in eachindex(x)]

    if typeof(K[1,1]) <: Number
        sK = Symmetric(K)
    else
        sK = K
    end
    return sK
end

"""
mean_log_marginal_loss(y_train, 𝒢::GP; add_constant=false)
# Description
Computes log marginal loss for each element in the output and averages the results.
Assumes noise-free observations.

log(p(y|X)) = -(1/2) * (y'*α + 2*sum(Diagonal(CK)) + n*log(2*pi))
where n is the number of training points and

# Arguments
- `y_train`: (Array). training outputs (prediction), must have the same number as x_train
- `𝒢`: (GP).
# Keyword Arguments
- `add_constant`: (bool). whether to give the exact value of the loss or leave out an added constant for efficiency.

"""
function mean_log_marginal_loss(y_train, 𝒢::GP; add_constant=false)
    n = length(𝒢.x_train)
    D = length(𝒢.x_train[1])

    ys = hcat(y_train...)' # n x D

    if add_constant
        c = sum([log(𝒢.CK[i,i]) for i in 1:n]) + 0.5*n*log(2*pi)
        total_loss=0.0
        for i in 1:D
            total_loss -= 0.5*ys[:,i]'*𝒢.α[:,i] + c
        end
    else
        total_loss=0.0
        for i in 1:D
            total_loss -= 0.5*ys[:,i]'*𝒢.α[:,i]
        end
    end

    return total_loss / D
end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection