text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

module OpticSim

import Unitful
using LinearAlgebra: eigen, svd, I, qr, dot, cross, norm, det, normalize, inv
using StaticArrays
using DataFrames: DataFrame
using Images
using Base: @.
using ForwardDiff
using StringEncodings

unzip(a) = map(x -> getfield.(a, x), fieldnames(eltype(a)))

include("GlassCat/GlassCat.jl")
import OpticSim.GlassCat: plot_indices, index, polyfit_indices, absairindex, absorption, info, glassid, glassname, glassforid, isair, findglass, modelglass, glassfromMIL, GlassID

include("constants.jl")
include("utilities.jl")
include("Geometry/Geometry.jl")
include("Optical/Optical.jl")
include("Visualization/Visualization.jl")
include("Examples/Examples.jl")
include("Optimization/Optimization.jl")

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

# This can be used to track NaN, particularly in ForwardDiff gradients, causing problems
# e.g. Diagnostics.testoptimization(lens = Examples.doubleconvex(NaNCheck{Float64}), samples = 1)

# struct NaNCheck{T<:Real} <: Real
#     val::T
#     function NaNCheck{T}(a::S) where {T<:Real, S<:Real}
#         @assert !(T <: NaNCheck)
#         new{T}(T(a))
#     end
# end
# export NaNCheck
# Base.isnan(a::NaNCheck{T}) where{T} = isnan(a.val)
# Base.isinf(a::NaNCheck{T}) where{T} = isinf(a.val)
# Base.typemin(::Type{NaNCheck{T}}) where{T} = NaNCheck{T}(typemin(T))
# Base.typemax(::Type{NaNCheck{T}}) where{T} = NaNCheck{T}(typemax(T))
# Base.eps(::Type{NaNCheck{T}}) where {T} = NaNCheck{T}(eps(T))
# Base.decompose(a::NaNCheck{T}) where {T} = Base.decompose(a.val)
# Base.round(a::NaNCheck{T}, m::RoundingMode) where {T} = NaNCheck{T}(round(a.val, m))

# struct NaNException <: Exception end

# # (::Type{Float64})(a::NaNCheck{S}) where {S<:Real} = NaNCheck{Float64}(Float64(a.val))
# (::Type{T})(a::NaNCheck{S}) where {T<:Integer,S<:Real} = T(a.val)
# (::Type{NaNCheck{T}})(a::NaNCheck{S}) where {T<:Real,S<:Real} = NaNCheck{T}(T(a.val))
# Base.promote_rule(::Type{NaNCheck{T}}, ::Type{T}) where {T<:Number} = NaNCheck{T}
# Base.promote_rule(::Type{T}, ::Type{NaNCheck{T}}) where {T<:Number} = NaNCheck{T}
# Base.promote_rule(::Type{S}, ::Type{NaNCheck{T}}) where {T<:Number, S<:Number} = NaNCheck{promote_type(T,S)}
# Base.promote_rule(::Type{NaNCheck{T}}, ::Type{S}) where {T<:Number, S<:Number} = NaNCheck{promote_type(T,S)}
# Base.promote_rule(::Type{NaNCheck{S}}, ::Type{NaNCheck{T}}) where {T<:Number, S<:Number} = NaNCheck{promote_type(T,S)}

# for op = (:sin, :cos, :tan, :log, :exp, :sqrt, :abs, :-, :atan, :acos, :asin, :log1p, :floor, :ceil, :float)
#     eval(quote
#         function Base.$op(a::NaNCheck{T}) where{T}
#             temp = NaNCheck{T}(Base.$op(a.val))
#             if isnan(temp)
#                 throw(NaNException())
#             end
#             return temp
#         end
#     end)
# end

# for op = (:+, :-, :/, :*, :^, :atan)
#     eval(quote
#         function Base.$op(a::NaNCheck{T}, b::NaNCheck{T}) where{T}
#             temp = NaNCheck{T}(Base.$op(a.val, b.val))
#             if isnan(temp)
#                 throw(NaNException())
#             end
#             return temp
#         end
#     end)
# end

# for op =  (:<, :>, :<=, :>=, :(==), :isequal)
#     eval(quote
#         function Base.$op(a::NaNCheck{T}, b::NaNCheck{T}) where{T}
#             temp = Base.$op(a.val, b.val)
#             return temp
#         end
#     end)
# end

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

end # module


text:
# This is code for computing a transit model and derivatives integrated over
# a time step, giving the fluence in units of time (since flux is normalized to unity).

# Now the version with derivatives:
function integrate_timestep_gradient!(param::Array{T,1},trans::Transit_Struct{T},t1::T,t2::T,tol::T,maxdepth::Int64,fint::Array{T,1}) where {T <: Real}
  # Keep track of maximum number of depths for recursion:
  depthmax = 0
  # Define third to save time in division:
  third = inv(3)
  # Keep track of how many evaluations of the transit function are carried out:
  neval = 0

  # When transit computation is carried out, we need to fill in the vector for the flux and derivatives
  # at the midpoint of the current refinement time range:
  function fill_flux!(tmid::T,fmid0::T,trans_mid::Transit_Struct{T},fmid::Array{T,1}) where {T <: Real}
#  fmid = Array{T}(undef,6+trans_mid.n)
  fmid[1] = fmid0   # flux
  fmid[2] = trans_mid.dfdrb[1]  # r derivative
  binv = inv(trans_mid.b)
  fmid[3] = trans_mid.dfdrb[2]*binv*param[2]^2*(param[1]-tmid)  # t0 derivative
  fmid[4] = trans_mid.dfdrb[2]*param[2]*binv*(tmid-param[1])^2  # v derivative
  fmid[5] = trans_mid.dfdrb[2]*param[3]*binv  # b0 derivative
  @inbounds for i=1:trans_mid.n+1
    fmid[5+i] = trans_mid.dfdg[i]  # d_i derivatives
  end
  return fmid
  end

  function fill_flux!(tmid::T,fmid0::T,trans_mid::Transit_Struct{T}) where {T <: Real}
  fmid = Array{T}(6+trans_mid.n)
  binv = inv(trans_mid.b)
  fmid[1] = fmid0   # flux
  fmid[2] = trans_mid.dfdrb[1]  # r derivative
  fmid[3] = trans_mid.dfdrb[2]*binv*param[2]^2*(param[1]-tmid)  # t0 derivative
  fmid[4] = trans_mid.dfdrb[2]*param[2]*binv*(tmid-param[1])^2  # v derivative
  fmid[5] = trans_mid.dfdrb[2]*param[3]*binv  # b0 derivative
  @inbounds for i=1:trans_mid.n+1
    fmid[5+i] = trans_mid.dfdg[i]  # d_i derivatives
  end
  return fmid
  end

  # This routine gives a simple model for the impact parameter versus time:
  function solver(time::T) where {T <: Real}
  # This predicts the impact parameter as a function of time in
  # terms of the transit_struct.
  # For now, transit is approximated as a straight line, where params are given as:
  # param[1] = t0 (mid-transit time; units of days/JD)
  # param[2] = v/R_* (sky velocity in terms of radius of the star - units are day^{-1}
  # param[3] = b (impact parameter at time t0)
  return sqrt(param[3]^2+(param[2]*(time-param[1]))^2)
  end

  # Function to define the vector integration in integrator:
#  function transit_flux_derivative!(tmid::T,fmid::Array{T,1}) where {T <: Real}
  function transit_flux_derivative(tmid::T) where {T <: Real}
#  if VERSION >=  v"0.7"
#    fmid = Array{T}(undef,6+trans.n)
#  else
#    fmid = Array{T}(6+trans.n)  # This line is taking a lot of time
#  end
  # Add the impact parameter to the transit structure:
  trans.b = solver(tmid)
  # Carry out the transit computation:
  fmid0 = transit_poly_g!(trans); neval += 1
  # Then, fillin the flux and derivatives in the fmid vector:
#  fill_flux!(tmid,fmid0,trans,fmid)  # So is this line
  fill_flux!(tmid,fmid0-1,trans,fmid)  # So is this line
#  println("tfd, r: ",trans.r," b: ",trans.b," f/df: ",fmid)
  return
  end

  # Function to define the vector integration in cubature:
  function transit_flux_derivative!(tmid::T,fmid::Array{T,1}) where {T <: Real}
  trans.b = solver(tmid)
  fmid0 = transit_poly_g!(trans); neval += 1
#  fmid .= fill_flux!(tmid,fmid0,trans)
  fill_flux!(tmid,fmid0-1,trans,fmid)
#  println("b: ",trans.b," f/df: ",fmid)
  return
  end

#  fint,ferr = hquadrature(trans.n+6,transit_flux_derivative,t1,t2,abstol=tol)
  # fint contains time-integrated flux (fluence) and derivatives:
# simpson(a::T, b::T, f::Function, I_of_f::Array{T,1}, i::T, epsilon::T, N::Int64, nf::Int64) where {T <: Real}
  tend = zero(T)
  fint .= simpson_vec(t1,t2,transit_flux_derivative!,fint,tend,tol,maxdepth,trans.n+6)
#  println("fint: ",fint)
#  fint2,ferr = hquadrature(trans.n+6,transit_flux_derivative,t1,t2,abstol=tol)
#  if maximum(abs,fint-fint2) > 1e-8
#    println("itg, r: ",trans.r," b: ",trans.b," f/df: ",fint," tol: ",tol," t2-t1: ",t2-t1)
#    println("cub, r: ",trans.r," b: ",trans.b," f/df: ",fint2)
#    println("f[t1]: ",transit_flux_derivative(t1)," f[t2]: ",transit_flux_derivative(t2))
#    read(STDIN,Char)
#  end
  # Return the number of evalutions and maximum depth for record-keeping:
return neval::Int64,depthmax::Int64
end

ext: jl


text:
"""
fractional flow formulation and solution of a Water injection process
for Corey-type oil and water relative permeabilities on a
Cartesian coordinate (not radial)

inputs [unit] \\(default):

   + muw [Pa.s] \\(1e-3)
   + muo [Pa.s] \\(2e-5)
   + ut [m/s] \\(1e-5)
   + phi [-] \\(0.2)
   + k [m^2] \\(1e-12)
   + swc [-] \\(0.1)
   + sor [-] \\(0.05)
   + kro0 [-] \\(0.9)
   + no [-] \\(1.8)
   + krw0 [-] \\(0.2)
   + nw [-] \\(4.2)
   + sw0 [-] \\(1.0)
   + sw_inj [-] \\(0.0)
   + L [m] \\(domain length = 1.0)
   + pv_inj [-] \\(injected pore volumes = 5)

Outputs:

   + xt_shock (x/t shock position)
   + sw_shock
   + xt_prf
   + sw_prf
   + t [s] \\time
   + p_inj [Pa] \\(injection pressure history)
   + R_oil [-] \\oil recovery factor

Usage:
   using Roots, Dierckx
   include("rel_perms.jl")

   + (xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R_oil)=frac_flow_wf()
   + (xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R_oil)=frac_flow_wf(muw, muo, ut, phi,
k, swc, sor, kro0, no, krw0, nw, sw0, sw_inj, L, pv_inj)
"""
function frac_flow_wf(;muw=1e-3, muo=2e-5, ut=1e-5, phi=0.2,
  k=1e-12, swc=0.1, sor=0.05, kro0=0.9, no=1.8, krw0=0.2,
  nw=4.2, sw0=0.0, sw_inj=1.0, L=1, pv_inj=5)
  # sws(sw::Real)=((sw>swc)*(sw<1-sor)*(sw-swc)/(1-sor-swc)+(sw>=1-sor)*1.0)
  # sws(sw::Array{Float64})=((sw.>swc).*(sw.<1-sor).*(sw-swc)/(1-sor-swc)+(sw.>=1-sor).*ones(size(sw)))
  # kro_new(sw)=((sw.>=swc).*kro0.*(1-sws(sw)).^no+(sw.<swc).*(1+(kro0-1)/swc*sw))
  # krw_new(sw)=((sw.<=1-sor).*krw0.*sws(sw).^nw+(sw.>1-sor).*(-(1-krw0)/sor.*(1.0-sw)+1.0))
  # dkrwdsw_new(sw)=(nw*krw0*(1/(1-sor-swc))*sws(sw).^(nw-1))
  # dkrodsw_new(sw)=((-kro0*no*(1-sws(sw)).^(no-1))/(-swc-sor+1))
  kro_new = sw -> kro(sw, kro0, sor, swc, no)
  krw_new = sw -> krw(sw, krw0, sor, swc, nw)
  dkrwdsw_new = sw -> dkrwdsw(sw, krw0, sor, swc, nw)
  dkrodsw_new = sw -> dkrodsw(sw, kro0, sor, swc, no)
  fw(sw)=((krw_new(sw)/muw)./(krw_new(sw)/muw+kro_new(sw)/muo))
  dfw(sw)=((dkrwdsw_new(sw)/muw.*(krw_new(sw)/muw+kro_new(sw)/muo)-
    (dkrwdsw_new(sw)/muw+dkrodsw_new(sw)/muo).*krw_new(sw)/muw)./
    (kro_new(sw)/muo+krw_new(sw)/muw).^2)
  eps1=1e-3
  dfw_num(sw) = (fw(sw+eps1)-fw(sw-eps1))/(2*eps1)
  ftot(sw)=kro_new(sw)/muo+krw_new(sw)/muw

# solve the nl equation to find the shock front saturation
  f_shock(sw)=(dfw(sw)-(fw(sw)-fw(sw0))/(sw-sw0))
  if f_shock(swc+eps())*f_shock(1-sor-eps())>0
      sw_shock = fzero(f_shock, (swc+1-sor)/2)
  else
      sw_shock = fzero(f_shock, [swc+eps(),1-sor-eps()])
  end
  s=collect(linspace(0.0,1.0,100))
  s1 = collect(linspace(sw_inj, sw_shock, 1000))
  xt_s1 = ut/phi*dfw.(s1)
  xt_s = ut/phi*dfw.(s)
  xt_shock = ut/phi*dfw(sw_shock)
  xt_prf=[xt_s1; xt_shock; xt_shock+eps(); 2*xt_shock]
  sw_prf=[s1; sw_shock; sw0; sw0]
  # println(xt_prf)
  # print(sw_shock)
# Go through the data first
  i=1
  while(true)
    if (i+1)==length(xt_prf)
      break
    elseif xt_prf[i]>=xt_prf[i+1]
      deleteat!(xt_prf, i+1)
      deleteat!(sw_prf, i+1)
    else
      i+=1
    end
  end
  # xt_prf_slim = unique(xt_prf)
  # sw_prf = sw_prf[indexin(xt_prf_slim, xt_prf)]
  # xt_prf = xt_prf_slim

# find the injection pressure history
  x = collect(linspace(0,L,1000))
  # sw_int = Spline1D([xt_prf; L/eps()], [sw_prf; sw0], k=1)
  sw_int = Spline1D(xt_prf, sw_prf, k=1)
  t_inj=pv_inj*phi*L/ut
  t = collect(linspace(0.0,t_inj, 200)) # [s] time
  p_inj = zeros(length(t))
  R_oil= zeros(length(t))
  p_inj[1]=trapz(x, ut./(k*(kro_new.(sw0*ones(size(x)))/muo+krw_new.(sw0*ones(size(x)))/muw)))
  for i in 2:length(t)
      xt = x/t[i]
      p_inj[i] = trapz(x, ut./(k*(kro_new.(sw_int(xt))/muo+krw_new.(sw_int(xt))/muw)))
      R_oil[i]=1.0-trapz(x/L, 1.0-sw_int(xt))/(1-sw0)
  end

# Return the results
    return xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R_oil
end

function trapz(x,y)
# copied from https://github.com/hwborchers/NumericalMath.jl/blob/master/src/integrate.jl
  n=length(x)
  if (length(y) != n)
    error("Vectors 'x', 'y' must be of same length")
  end
  r=0
  for i in 2:n
    r+=(x[i] - x[i-1])*(y[i] + y[i-1])
  end
  return r/2.0
end

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
using Soss
using Test
using MeasureTheory
using TransformVariables
using Aqua
Aqua.test_all(Soss; ambiguities=false, unbound_args=false)

include("examples-list.jl")

@testset "Soss.jl" begin
    @testset "Unit tests" begin
        @testset "Linear model" begin
            # include("linear-model.jl")
        end

        @testset "Transforms" begin
            include("transforms.jl")
        end
    end

    @testset "Examples" begin
        for example in EXAMPLES
            @testset "Run example: $(example[1])" begin
                # example_file = joinpath(EXAMPLESROOT, "example-$(example[2]).jl")
                # extra_example_tests = joinpath(TESTROOT, "extra-example-tests", "$(example[2]).jl")
                # @info("Running $(example_file)")
                # include(example_file)
                # if isfile(extra_example_tests)
                #     @info("Running $(extra_example_tests)")
                #     include(extra_example_tests)
                # end
            end
        end
    end

    @testset "Nested models" begin
        inner = @model a, b begin
            p ~ Beta(a, b)
            x ~ Normal(p, 1.0) |> iid(3)
            return x
        end

        outer = @model sub begin
            a ~ Beta(0.5, 0.5)
            b ~ Beta(1, 0.5)
            m ~ sub(a = a, b = b)
        end

        x = rand(outer(sub=inner)).m
        post = outer(sub=inner) | (m = (x=x,),)
        t = xform(post)
        @test logdensity(post, transform(t, randn(3))) isa Float64
    end

    @testset "Predict" begin
        m = @model begin
            p ~ Uniform()
            y ~ Bernoulli(p)
            return y
        end

        mean(predict(m(), [(p=p,) for p in rand(10000)])) isa Float64
    end

    @testset "https://github.com/cscherrer/Soss.jl/issues/258" begin
        m1 = @model begin
            x1 ~ Soss.Normal(0.0, 1.0)
            x2 ~ Dists.LogNormal(0.0, 1.0)
            return x1^2/x2
        end

        m2 = @model m begin
            μ ~ m
            y ~ Soss.Normal(μ, 1.0)
        end

        mm = m2(m=m1())

        @test xform(mm|(y=1.0,)) isa TransformVariables.TransformTuple
        @test basemeasure(mm | (y=1.0,)) isa ProductMeasure
        @test testvalue(mm) isa NamedTuple
    end

    @testset "https://github.com/cscherrer/Soss.jl/issues/258#issuecomment-819035325" begin
        m1 = @model begin
            x1 ~ Soss.Normal(0.0, 1.0)
            x2 ~ Dists.MvNormal(fill(x1,2), ones(2))
            return x2
        end

        m2 = @model m begin
            μ ~ m
            y ~ For(μ) do x
                Soss.Normal(x, 1.0)
            end
        end

        mm = m2(m=m1())

        @test xform(mm|(;y=1.0,)) isa TransformVariables.TransformTuple
        @test basemeasure(mm | (y=1.0,)) isa ProductMeasure
        @test testvalue(mm) isa NamedTuple
    end

    @testset "https://github.com/cscherrer/Soss.jl/issues/305" begin
        m = @model begin
            x ~ For(3) do j Normal(μ=j) end
        end;

        @test logpdf(m(), rand(m())) isa Float64
    end

    @testset "Local variables" begin
        # https://github.com/cscherrer/Soss.jl/issues/253

        m = @model begin
            a ~ For(3) do x Normal(μ=x) end
            x ~ Normal(μ=sum(a))
        end

        @test digraph(m).N == Dict(:a => Set([:x]), :x => Set())
    end

    @testset "Doctests" begin
        include("doctests.jl")
    end

    @testset "Distributions" begin
        m = @model begin
            a ~ Normal() |> iid(3)
            b ~ Dists.Normal() |> iid(3)
            c ~ For(3) do i
                Normal(μ = a[i] +b[i])
            end
        end

        c = rand(m()).c

        post = m() | (c=c,)

        @test transform(xform(post), randn(6)) isa NamedTuple

        @testset "logdensity" begin
            dat = randn(100)
            m = Soss.@model n begin
                μ ~ Dists.Normal()
                σ ~ Dists.Exponential()
                data ~ Dists.Normal(μ, σ) |> iid(n)
                return (; data)
            end
            mod = m( (; n = length(dat) ) )
            post = mod | (data = dat,)

            @test logdensity( mod( (μ = 1., σ = 2., data = dat) ) ) == logdensity( post( (μ = 1., σ = 2.) ) )
        end

    end

end


text:
@testset "CURAND" begin

if !isdefined(CuArrays, :CURAND)
@warn "Not testing CURAND"
else
using CuArrays.CURAND
@info "Testing CURAND $(CURAND.version())"

CURAND.seed!()

# in-place
for (f,T) in ((rand!,Float32),
              (randn!,Float32),
              (rand_logn!,Float32),
              (rand_poisson!,Cuint)),
    d in (2, (2,2), (2,2,2))
    A = CuArray{T}(undef, d)
    f(A)
end

# out-of-place, with implicit type
for (f,T) in ((curand,Float32), (curandn,Float32), (curand_logn,Float32),
              (curand_poisson,Cuint),(rand,Float64), (randn,Float64)),
    args in ((2,), (2, 2))
    A = f(args...)
    @test eltype(A) == T
end

# out-of-place, with type specified
for (f,T) in ((curand,Float32), (curandn,Float32), (curand_logn,Float32),
              (curand,Float64), (curandn,Float64), (curand_logn,Float64),
              (curand_poisson,Cuint), (rand,Float32), (randn,Float32),
              (rand_logn,Float32), (rand,Float64), (randn,Float64),
              (rand_logn,Float64), (rand_poisson,Cuint)),
    args in ((T, 2), (T, 2, 2), (T, (2, 2)))
    A = f(args...)
    @test eltype(A) == T
end

# unsupported types that fall back to GPUArrays
for (f,T) in ((curand,Int64),),
    args in ((T, 2), (T, 2, 2), (T, (2, 2)))
    A = f(args...)
    @test eltype(A) == T
end
for (f,T) in ((rand!,Int64),),
    d in (2, (2,2), (2,2,2))
    A = CuArray{T}(undef, d)
    f(A)
end

@test_throws ErrorException randn!(CuArray{Cuint}(undef, 10))
@test_throws ErrorException rand_logn!(CuArray{Cuint}(undef, 10))
@test_throws ErrorException rand_poisson!(CuArray{Float64}(undef, 10))

end

end


text:
using GLPlot, GLVisualize, GLAbstraction, Colors, GeometryTypes, FileIO
using Reactive, GLWindow, NIfTI
const If0 = GLVisualize.Intensity{1, Float32}
window = GLPlot.init()

layout = [
    ["Volume Plot", "Slice X"],
    ["Slice Y", "Slice Z"]
]

screens = GLVisualize.layoutscreens(window, layout)
function volume_data(N)
    vol = Float32[sin(x/15.0) + sin(y/15.0) + sin(z/15.0) for x=1:N, y=1:N, z=1:N]
    max = maximum(vol); min     = minimum(vol)
    vol = (vol .- min) ./ (max .- min)
end

vol = volume_data(128)

# to change color_map use the interactive edit menu or pass an array of color
# via color_map = Vector{Colorant}
p1 = glplot(vol, screen = screens[2][1])
slice_screens = vcat(screens[2][2], screens[1]...)
for i = 1:3
    range_s = play_widget(1:size(vol, i))
    slice = map(range_s) do slice_idx
        idx = ntuple(d-> d == i ? slice_idx : (:), Val{3})
        # This conversion is necessary but will be automatic soon!
        If0.(view(vol, idx...))
    end
    plot = glplot(
        slice, screen = slice_screens[i],
        color_norm = Vec2f0(0, 1),
        camera = :orthographic_pixel
    )
    center!(slice_screens[i], :orthographic_pixel)
end


text:
__precompile__()
"""
Module with reachability algorithms.
"""
module ReachSets

using ..Utils
using LazySets, MathematicalSystems, MathematicalPredicates, HybridSystems,
      Expokit, Optim, LinearAlgebra, SparseArrays
import IntervalMatrices
using IntervalMatrices: IntervalMatrix

import ProgressMeter
using ProgressMeter: Progress
# no-op progress meter for unbounded time
function ProgressMeter.update!(::Nothing, ::Int) end

# fix namespace conflicts
using LazySets: LinearMap, AffineMap, ResetMap, Interval
using LazySets.Approximations: AbstractDirections

import LazySets: set
import LazySets.Approximations: project

using Reachability: info, warn,
                    @timing,
                    Options, OptionSpec, TwoLayerOptions, AbstractOptions,
                    validate_and_wrap_options, print_option_spec,
                    haskey_specified,
                    validate_solver_options_and_add_default_values!

# ========================================
# Discretize and compute bloating factors
# ========================================
include("discretize.jl")

# ==========================
# Property checking results
# ==========================
include("CheckSolution.jl")

export CheckSolution

# =============================
# Property checking algorithms
# =============================

# dictionary of registered algorithms
available_algorithms_check = Dict{String, Dict{String, Any}}()

include("ContinuousPost/BFFPSV18/check_blocks.jl")
include("ContinuousPost/BFFPSV18/check_property.jl")
include("ContinuousPost/BFFPSV18/partitions.jl")
include("ContinuousPost/BFFPSV18/compute_dimensions.jl")

# "explicit" backends
push!(available_algorithms_check, "explicit_blocks"=>Dict("func"=>check_blocks,
                                                          "is_explicit"=>true))

export available_algorithms_check,
       check_property

# ====================================================
# Algorithms to find a threshold for property checking
# ====================================================
include("tune.jl")
export tune_δ

# =====================
# Reachability results
# =====================
include("ReachSet.jl")
include("Flowpipe.jl")
include("ReachSolution.jl")

export AbstractReachSet, ReachSet, SparseReachSet,
       set, time_start, time_end,
       Flowpipe,
       ReachSolution,
       n_sets

# ===============
# Post operators
# ===============
include("PostOperator.jl")
include("ContinuousPost/ContinuousPost.jl")
include("DiscretePost/DiscretePost.jl")

export PostOperator,
       ContinuousPost,
       DiscretePost,
       init, init!,
       post,
       tube⋂inv

# ==========================
# Continuous post operators
# ==========================
include("ContinuousPost/BFFPSV18/BFFPSV18.jl")
include("ContinuousPost/BFFPSV18/reach.jl")
include("ContinuousPost/BFFPSV18/reach_blocks.jl")
include("ContinuousPost/BFFPSV18/reach_blocks_wrapping_effect.jl")

include("ContinuousPost/GLGM06/GLGM06.jl")

include("ContinuousPost/TMJets/TMJets.jl")

include("ContinuousPost/BFFPS19/BFFPS19.jl")

include("ContinuousPost/ASB07/ASB07.jl")

include("ContinuousPost/ASB07_decomposed/ASB07_decomposed.jl")

# ========================
# Reachability Algorithms
# ========================
import Reachability.check_aliases_and_add_default_value!

# dictionary of registered algorithms
available_algorithms = Dict{String, Dict{String, Any}}()

# "explicit" backends
push!(available_algorithms, "explicit_blocks"=>Dict("func"=>reach_blocks!,
                                                    "is_explicit"=>true))

push!(available_algorithms,
      "wrap"=>Dict("func"=>reach_blocks_wrapping_effect!,
                   "is_explicit"=>true))

export available_algorithms

# ==========================
# Discrete post operators
# ==========================
include("DiscretePost/LazyDiscretePost.jl")
include("DiscretePost/ConcreteDiscretePost.jl")
include("DiscretePost/DecomposedDiscretePost.jl")

# ==============================================
# Projection of the reach set in two dimensions
# ==============================================
include("project_reach.jl")

export project_reach

end # module ReachSets


text:
push!(LOAD_PATH, ".")
macro tnum()
	return :(testnum += 1)
end
using Lexer
using Error
using HPInt
using Cairo
using Images
testnum = 0
function testNum(num)
	return string(num) * ". "
end
function parseT(str)
	HPInt.parse(Lexer.lex(str))
end
function analyze(str)
	HPInt.analyze(parseT(str))
end
function interpret(str)
	HPInt.calc(analyze(str))
end

test1="(with ((x 1)) x)"
test2="((lambda (x) x) 1)"
println(testNum(@tnum()) * "analyze " * test1)
try
	HPInt.analyze(HPInt.parse(Lexer.lex(test1)))
catch
	println("Failed test " * testNum(testnum))
	println(HPInt.analyze(HPInt.parse(Lexer.lex(test1))))
end
println(testNum(@tnum()) * "analyze " * test2)
try
	HPInt.analyze(HPInt.parse(Lexer.lex(test2)))
catch
	println("Failed test " * testNum(testnum))
	println(HPInt.analyze(HPInt.parse(Lexer.lex(test2))))
end
try
	if(string(HPInt.analyze(HPInt.parse(Lexer.lex(test1))))==string(analyze(HPInt.parse(Lexer.lex(test2)))))
		println("Passed tests " * testNum(testnum-1) * "and " * testNum(testnum))
	else
		println("Failed test " * testNum(testnum-1) * "and " * testNum(testnum) * "These should be the same:")
		println(HPInt.analyze(HPInt.parse(Lexer.lex(test1))))
		println(HPInt.analyze(HPInt.parse(Lexer.lex(test2))))
	end
catch
	println("Failed test " * testNum(testnum-1) * "and " * testNum(testnum) * "These should be the same:")
	println(HPInt.analyze(HPInt.parse(Lexer.lex(test1))))
	println(HPInt.analyze(HPInt.parse(Lexer.lex(test2))))
end
println(testNum(@tnum()) * "calc " * test1)
try
	HPInt.calc(test1)
catch
	println("Failed test " * testNum(testnum))
	println(HPInt.calc(test1))
end
println(testNum(@tnum()) * "calc " * test2)
try
	HPInt.calc(test2)
catch
	println("Failed test " * testNum(testnum))
	println(HPInt.calc(test2))
end
try
	if(string(HPInt.calc(test1))==string(HPInt.calc(test2)))
		println("Passed tests " * testNum(testnum-1) * "and " * testNum(testnum))
	else
		println("Failed test " * testNum(testnum-1) * "and " * testNum(testnum) * "These should be the same:")
		println(HPInt.calc(test1))
		println(HPInt.calc(test2))
	end
catch
	println("Failed test " * testNum(testnum-1) * "and " * testNum(testnum) * "These should be the same:")
	println(HPInt.calc(test1))
	println(HPInt.calc(test2))
end

test1="(with ((x 1) (y 2)) y)"
test2="(with ((x 1) (y 2)) y)"
println(testNum(@tnum()) * "analyze " * test1)
try
	HPInt.analyze(HPInt.parse(Lexer.lex(test1)))
catch
	println("Failed test " * testNum(testnum))
	println(HPInt.analyze(HPInt.parse(Lexer.lex(test1))))
end
println(testNum(@tnum()) * "analyze " * test2)
try
	HPInt.analyze(HPInt.parse(Lexer.lex(test2)))
catch
	println("Failed test " * testNum(testnum))
	println(HPInt.analyze(HPInt.parse(Lexer.lex(test2))))
end
try
	if(string(HPInt.analyze(HPInt.parse(Lexer.lex(test1))))==string(HPInt.analyze(HPInt.parse(Lexer.lex(test2)))))
		println("Passed tests " * testNum(testnum-1) * "and " * testNum(testnum))
	else
		println("Failed test " * testNum(testnum-1) * "and " * testNum(testnum) * "These should be the same:")
		println(HPInt.analyze(HPInt.parse(Lexer.lex(test1))))
		println(HPInt.analyze(HPInt.parse(Lexer.lex(test2))))
	end
catch
	println("Failed test " * testNum(testnum-1) * "and " * testNum(testnum) * "These should be the same:")
	println(HPInt.analyze(HPInt.parse(Lexer.lex(test1))))
	println(HPInt.analyze(HPInt.parse(Lexer.lex(test2))))
end
println(testNum(@tnum()) * "calc " * test1)
try
	HPInt.calc(test1)
catch
	println("Failed test " * testNum(testnum))
	println(HPInt.calc(test1))
end
println(testNum(@tnum()) * "calc " * test2)
try
	HPInt.calc(test2)
catch
	println("Failed test " * testNum(testnum))
	println(HPInt.calc(test2))
end
try
	if(string(HPInt.calc(test1))==string(HPInt.calc(test2)))
		println("Passed tests " * testNum(testnum-1) * "and " * testNum(testnum))
	else
		println("Failed test " * testNum(testnum-1) * "and " * testNum(testnum) * "These should be the same:")
		println(HPInt.calc(test1))
		println(HPInt.calc(test2))
	end
catch
	println("Failed test " * testNum(testnum-1) * "and " * testNum(testnum) * "These should be the same:")
	println(HPInt.calc(test1))
	println(HPInt.calc(test2))
end

println(testNum(@tnum()) * "Load and Save")
test="(with ((cat (simple_load \"cat_256.png\"))) (simple_save cat \"output"*string(testnum)*".png\"))"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create image from sample code.")
    	interpret(test)
end
println(testNum(@tnum()) * "Emboss")
test="(with ((cat (simple_load \"cat_256.png\"))) (simple_save (emboss cat) \"output"*string(testnum)*".png\"))"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create image from sample code.")
    	interpret(test)
end
println(testNum(@tnum()) * "Inner Shadow")
test="(with ((cat (simple_load \"cat_256.png\"))) (simple_save (inner_shadow cat) \"output"*string(testnum)*".png\"))"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create image from sample code.")
    	interpret(test)
end
println(testNum(@tnum()) * "Drop Shadow")
test="(with ((cat (simple_load \"cat_256.png\"))) (simple_save (drop_shadow cat) \"output"*string(testnum)*".png\"))"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create image from sample code.")
    	interpret(test)
end
println(testNum(@tnum()) * "Max")
test="(with ((cat (simple_load \"cat_256.png\"))) (simple_save (max 0.5 cat) \"output"*string(testnum)*".png\"))"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create image from sample code.")
    	interpret(test)
end
println(testNum(@tnum()) * "Min")
test="(with ((cat (simple_load \"cat_256.png\"))) (simple_save (min 0.5 cat) \"output"*string(testnum)*".png\"))"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create image from sample code.")
    	interpret( test)
end
println(testNum(@tnum()) * "Decrement")
test="(with ((cat (simple_load \"cat_256.png\"))) (simple_save (- 1 cat) \"output"*string(testnum)*".png\"))"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create image from sample code.")
    	interpret(test)
end
println(testNum(@tnum()) * "Addition")
test="(with ((cat (simple_load \"cat_256.png\"))) (simple_save (+ cat cat) \"output"*string(testnum)*".png\"))"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create image from sample code.")
    	interpret(test)
end
println(testNum(@tnum()) * "Subtraction")
test="(with ((cat (simple_load \"cat_256.png\"))) (simple_save (- cat cat) \"output"*string(testnum)*".png\"))"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create image from sample code.")
    	interpret( test)
end
println(testNum(@tnum()) * "Render Text")
test="(with ((cat (simple_load \"cat_256.png\"))) (simple_save (min (render_text \"CAT\" 40 230) cat) \"output"*string(testnum)*".png\"))"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create image from sample code.")
    	interpret(test)
end
println(testNum(@tnum()) * "Complex")
test="(with ((cat (simple_load \"cat_256.png\"))) (simple_save (* (- 1 cat) cat) \"output"*string(testnum)*".png\"))"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create image from sample code.")
    	interpret(test)
end

println(testNum(@tnum()) * "Swirl Load")
test="(with ((cat (simple_load \"swirl_256.png\"))) (simple_save (* (- 1 cat) cat) \"output"*string(testnum)*".png\"))"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create image from sample code.")
    	interpret(test)
end
println(testNum(@tnum()) * "Sample Hello World")
test="(with ( (base_img (render_text \"Hello\" 25 100)) (swirl (simple_load \"swirl_256.png\")) ) (with ( (ds (drop_shadow base_img)) ) (with ( (tmp4 (+ (* (+ (min ds base_img) (- 1 base_img)) base_img) (* (- 1 base_img) swirl) )) ) (with ( (tmp5 (- 1 (emboss tmp4))) (base_img2 (render_text \"world!\" 5 200)) ) (with ( (is (inner_shadow base_img2)) ) (with ( (tmp6 (max base_img2 (* (- 1 base_img2) is))) ) (with ( (output (min tmp5 tmp6)) ) (simple_save output \"output"*string(testnum)*".png\") ) ) ) ) ) ) )"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create image from sample code.")
    	interpret(test)
end

println(testNum(@tnum()) * "Cat picture")
test="(with ((cat (simple_load \"cat_256.png\"))) (with ((inv (- 1 cat))) (with ( (output (drop_shadow inv)) ) (simple_save output \"output"*string(testnum)*".png\"))))"
try
	interpret(test)
	println("Passed test " * testNum(testnum))
catch
	println("Failed to create cat image.")
	interpret(test)
end

max_stars_repo_name: scottcrossen/JuliaDevelopment

max_issues_repo_name: scottcrossen/JuliaDevelopment

max_forks_repo_name: scottcrossen/JuliaDevelopment


text:
"""
    module ArnoldiMethodTransformations

Provides convenience wrapper for interfacing with the package ArnoldiMethod.

Implements the shift-and-invert transformation detailed [here](https://haampie.github.io/ArnoldiMethod.jl/stable/).

The main functions are `partialschur(A,[B],σ; kwargs...)` and `partialeigen(A,σ; kwargs...)`

The constants `USOLVER`, `PSOLVER`, and `MSOLVER` are exported.
"""
module ArnoldiMethodTransformations

export PSOLVER
export PSolver
export MSOLVER
export MSolver
export USOLVER
export USolver

using ArnoldiMethod
using LinearAlgebra
using LinearMaps
#using MUMPS
using Pardiso
using SparseArrays

abstract type AbstractSolver end
struct PSolver <: AbstractSolver end
struct MSolver <: AbstractSolver end
struct USolver <: AbstractSolver end

# solver constants (exported)
const PSOLVER = PSolver()
const MSOLVER = MSolver()
const USOLVER = USolver()

# colors for pretty printing
const DEFAULT_SOLVER = USOLVER
const SOLVER_COLOR = :cyan
const PRINTED_TYPE_COLOR = 171

# main structure used internally
struct ShiftAndInvert{M,TLU,TB,V,Σ,TZ}
    LU::TLU
    B::TB
    temp1::V
    temp2::V
    σ::Σ
    diag_inv_B::Bool
    issymmetric::Bool
    Z::TZ
end

# for ordinary EVP with Dense Matrix
function ShiftAndInvert(
            A::Matrix{TA};
            σ::Σ=0,
            lupack::AbstractSolver = DEFAULT_SOLVER
            ) where {TA,TB,Σ<:Number}

    diag_inv_B=true
    T = promote_type(TA,Σ)
    α = A-σ*I
    β = Matrix(one(T)I,size(A))

    temp1 = Vector{T}(undef,size(α,1))
    temp2 = Vector{T}(undef,size(α,1))
    issym = issymmetric(α) & issymmetric(β)
    Z = spzeros(T,size(β)...)

    M, a_lu = initialize_according_to_package(lupack,α,issym,T,temp1,temp2)
    return ShiftAndInvert{typeof(M),typeof(a_lu),typeof(β),typeof(temp1),typeof(σ),typeof(Z)}(a_lu,β,temp1,temp2,σ,diag_inv_B,issym,Z)
end
# for generalized EVP with Dense Matrices
function ShiftAndInvert(
            A::Matrix{TA},
            B::Matrix{TB};
            σ::Σ=0,
            diag_inv_B::Bool=false,
            lupack::AbstractSolver = DEFAULT_SOLVER
            ) where {TA,TB,Σ<:Number}

    T = promote_type(TA,TB,Σ)
    if diag_inv_B
        for i ∈ diagind(B) B[i]=1/B[i] end
        α = B*A-σ*I
        for i ∈ diagind(B) B[i]=1/B[i] end
        β = Matrix(one(T)I,size(B))
    else
        α = A-σ*B
        β = convert(Matrix{T},B)
    end
    temp1 = Vector{T}(undef,size(α,1))
    temp2 = Vector{T}(undef,size(α,1))
    issym = issymmetric(α) & issymmetric(β)
    Z = spzeros(T,size(β)...)

    M, a_lu = initialize_according_to_package(lupack,α,issym,T,temp1,temp2)
    return ShiftAndInvert{typeof(M),typeof(a_lu),typeof(β),typeof(temp1),typeof(σ),typeof(Z)}(a_lu,β,temp1,temp2,σ,diag_inv_B,issym,Z)
end
# for ordinary EVP with Sparse Matrix
function ShiftAndInvert(
            A::SparseMatrixCSC{TA};
            σ::Σ=0,
            lupack::AbstractSolver = DEFAULT_SOLVER
            ) where {TA,TB,Σ}

    diag_inv_B=true
    T = promote_type(TA,Σ)
    α = A-σ*I
    β = sparse(one(T)I,size(A))
    temp1 = Vector{T}(undef,size(α,1))
    temp2 = Vector{T}(undef,size(α,1))
    issym = issymmetric(α) & issymmetric(β)
    Z = spzeros(T,size(β)...)

    # initialize according to package used
    M, a_lu = initialize_according_to_package(lupack,α,issym,T,temp1,temp2)
    return ShiftAndInvert{typeof(M),typeof(a_lu),typeof(β),typeof(temp1),typeof(σ),typeof(Z)}(a_lu,β,temp1,temp2,σ,diag_inv_B,issym,Z)
end
# for generalized EVP with Sparse Matrices
function ShiftAndInvert(
            A::SparseMatrixCSC{TA},
            B::SparseMatrixCSC{TB};
            σ::Σ=0,
            diag_inv_B::Bool=false,
            lupack::AbstractSolver = DEFAULT_SOLVER
            ) where {TA,TB,Σ}

    T = promote_type(TA,TB,Σ)
    if diag_inv_B
        for i ∈ diagind(B) B[i]=1/B[i] end
        α = B*A-σ*I
        for i ∈ diagind(B) B[i]=1/B[i] end
        β = sparse(one(T)I,size(B))
    else
        α = A-σ*B
        β = convert(SparseMatrixCSC{T,Int},B)
    end
    temp1 = Vector{T}(undef,size(α,1))
    temp2 = Vector{T}(undef,size(α,1))
    issym = issymmetric(α) & issymmetric(β)
    Z = spzeros(T,size(β)...)

    # initialize according to package used
    M, a_lu = initialize_according_to_package(lupack,α,issym,T,temp1,temp2)
    return ShiftAndInvert{typeof(M),typeof(a_lu),typeof(β),typeof(temp1),typeof(σ),typeof(Z)}(a_lu,β,temp1,temp2,σ,diag_inv_B,issym,Z)
end
# if not both dense or both sparse
ShiftAndInvert(A,B;kwargs...) = throw("typeof(A)=$(typeof(A)), typeof(B)=$(typeof(B)). Either both sparse or both dense")

# define action of ShiftAndInvert
function (SI::ShiftAndInvert{M})(y,x) where M<:Union{USolver,MSolver}
    #first multiply Bv
    if SI.diag_inv_B
        for i ∈ eachindex(SI.temp1) SI.temp1[i] = x[i] end
    else
        mul!(SI.temp1, SI.B, x)
    end
    # then left-divide (A-σB)⁻¹*(Bv)
    ldiv!(y, SI.LU, SI.temp1)
    return nothing
end
function (SI::ShiftAndInvert{PSolver})(y,x)
    #first multiply Bv
    if SI.diag_inv_B
        for i ∈ eachindex(SI.temp1) SI.temp1[i] = x[i] end
    else
        mul!(SI.temp1, SI.B, x)
    end
    # then left-divide (A-σB)⁻¹*(Bv)
    pardiso(SI.LU,SI.temp2,SI.Z,SI.temp1)
    for i ∈ eachindex(y) y[i] = SI.temp2[i] end
    return nothing
end

# generate a solver-dependent LU object
# UMFPACK
initialize_according_to_package(::USolver,α,args...) = USOLVER, lu(α)
# MUMPS
function initialize_according_to_package(::MSolver,α,args...)
    throw("MUMPS still under construction")
    try MPI catch; throw(ErrorException("MUMPS not loaded. Try again after `using MUMPS`")) end
    MPI.Initialized() ? nothing : MPI.Init()
    return MSOLVER, mumps_factorize(α)
end
# PARDISO
function initialize_according_to_package(::PSolver,α,issym,type,x,y)
    try PardisoSolver catch; throw(ErrorException("Pardiso not loaded. Try again after `using Pardiso`")) end
    A_lu = PardisoSolver()
    set_iparm!(A_lu,1,1) # don't revert to defaults
    set_iparm!(A_lu,12,1) # transpose b/c of CSR vs SCS
    # x = Vector{type}(undef,size(α,1))
    # y = Vector{type}(undef,size(α,1))
    if issym & (type<:Real)
        set_matrixtype!(A_lu,Pardiso.REAL_SYM)
        pardiso(A_lu,x,sparse(tril(α)),y)
    elseif issym & (type<:Complex)
        set_matrixtype!(A_lu,Pardiso.COMPLEX_SYM)
        pardiso(A_lu,x,sparse(tril(α)),y)
    elseif !issym & (type<:Real)
        set_matrixtype!(A_lu,Pardiso.REAL_NONSYM)
        pardiso(A_lu,x,sparse(α),y)
    else # if !issym & type<:Complex
        set_matrixtype!(A_lu,Pardiso.COMPLEX_NONSYM)
        pardiso(A_lu,x,sparse(α),y)
    end
    set_phase!(A_lu,12) # analyze and factorize
    set_phase!(A_lu,33) # set to solve for future calls
    return PSOLVER, A_lu
end

# pretty printing
function Base.show(io::IO,si::ShiftAndInvert{M}) where M
    printstyled(io,"ShiftAndInvert ",color=PRINTED_TYPE_COLOR)
    print(io,"(using ")
    if M==MSolver
        printstyled(io,"MUMPS",color=SOLVER_COLOR)
    elseif M==PSolver
        printstyled(io,"Pardiso",color=SOLVER_COLOR)
    elseif M==USolver
        printstyled(io,"UMFPACK",color=SOLVER_COLOR)
    end
    print(io,")")
end

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

# add method to partialschur(::ShiftAndInvert)
function ArnoldiMethod.partialschur(si::ShiftAndInvert; kwargs...)
    a = LinearMap{eltype(si.B)}(si, size(si.B,1); ismutating=true, issymmetric=si.issymmetric)
    return partialschur(a; kwargs...)
end
# add method to partialschur(A,σ)
function ArnoldiMethod.partialschur(
            A::AbstractArray,
            σ::Number;
            lupack::AbstractSolver=DEFAULT_SOLVER,
            kwargs...)
    return partialschur(ShiftAndInvert(A; σ=σ, lupack=lupack); kwargs...)
end
# add method to partialschur(A,B,σ)
function ArnoldiMethod.partialschur(
            A::AbstractArray,
            B::AbstractArray,
            σ::Number;
            diag_inv_B::Bool = isdiag(B) && !any(iszero.(diag(B))),
            lupack::AbstractSolver=DEFAULT_SOLVER,
            kwargs...)
    return partialschur(ShiftAndInvert(A, B; σ=σ, diag_inv_B=diag_inv_B, lupack=lupack); kwargs...)
end

# add method to partialeigen(decomp,σ)
"""
    partialeigen(decomp, σ) -> λ, v

    Transforms a partial Schur decomposition into an eigendecomposition, outputting evals `λ` and evecs `v`.
    It undoes the shift-and-invert of the eigenvalues by `σ`.
"""
function ArnoldiMethod.partialeigen(
            decomp::ArnoldiMethod.PartialSchur,
            σ::Number)
    λ, v = partialeigen(decomp)
    foreach(i->λ[i]=σ+1/λ[i],eachindex(λ))
    return λ, v
end

################################################################################
"""
    struct ShiftAndInvert

Container for arrays in `ArnoldiMethod.partialschur`

-------------

    function ShiftAndInvert(A, [B]; σ=0, diag_inv_B=false, lupack=USOLVER) -> si

create a LinearMap object to feed to ArnoldiMethod.partialschur which transforms `Ax=λBx` into `(A-σB)⁻¹Bx=x/(λ-σ)`.

Set `diag_inv_B=true` if `B` is both diagonal and invertible, so that it is easy to compute `B⁻¹`. In this case instead return a linear map `(B⁻¹A-σI)⁻¹`, which has same evals as above.

`A` and `B` must both be sparse or both dense. `A`, `B`, `σ` need not have common element type.

Keyword `lupack` determines what linear algebra library to use. Options are `PSOLVER`, `MSOLVER`, `USOLVER` (the default).
For all solvers besides UMFPACK, the appropriate library must be loaded at the top level.
----------

    function ::ShiftAndInvert(y,x)

`A\\B*x = y`
"""
ShiftAndInvert

"""
    partialschur(A, [B], σ; [diag_inv_B=false, lupack=USOLVER, kwargs...]) -> decomp, history

Partial Schur decomposition of `A`, with shift `σ` and mass matrix `B`, solving `Av=σBv`

Keyword `diag_inv_B` defaults to `true` if `B` is both diagonal and invertible. This enables
a simplified shift-and-invert scheme.

Keyword `lupack` determines what linear algebra library to use. Options are `PSOLVER`, `MSOLVER`, `USOLVER` (the default).
For all solvers besides UMFPACK, the appropriate library must be loaded at the top level.

For other keywords, see ArnoldiMethod.partialschur

see also: [`partialeigen`](@ref) in ArnoldiMethod
"""
partialschur

end # module


text:
module Bridges

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

const MOIB = MOI.Bridges # used in macro

const SVF = MOI.SingleVariable
const VVF = MOI.VectorOfVariables
const SAF{T} = MOI.ScalarAffineFunction{T}
const VAF{T} = MOI.VectorAffineFunction{T}
const SQF{T} = MOI.ScalarQuadraticFunction{T}
const VQF{T} = MOI.VectorQuadraticFunction{T}

const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

include("bridge.jl")
include("bridgeoptimizer.jl")
include("singlebridgeoptimizer.jl")
include("lazybridgeoptimizer.jl")

"""
    full_bridge_optimizer(model::MOI.ModelLike, ::Type{T}) where T

Returns a `LazyBridgeOptimizer` bridging `model` for every bridge defined in
this package and for the coefficient type `T`.
"""
function full_bridge_optimizer(model::MOI.ModelLike, ::Type{T}) where T
    bridged_model = LazyBridgeOptimizer(model)
    add_bridge(bridged_model, GreaterToLessBridge{T})
    add_bridge(bridged_model, LessToGreaterBridge{T})
    add_bridge(bridged_model, NonnegToNonposBridge{T})
    add_bridge(bridged_model, NonposToNonnegBridge{T})
    add_bridge(bridged_model, ScalarizeBridge{T})
    add_bridge(bridged_model, VectorizeBridge{T})
    add_bridge(bridged_model, ScalarSlackBridge{T})
    add_bridge(bridged_model, VectorSlackBridge{T})
    add_bridge(bridged_model, ScalarFunctionizeBridge{T})
    add_bridge(bridged_model, VectorFunctionizeBridge{T})
    add_bridge(bridged_model, SplitIntervalBridge{T})
    add_bridge(bridged_model, QuadtoSOCBridge{T})
    add_bridge(bridged_model, GeoMeanBridge{T})
    add_bridge(bridged_model, SquareBridge{T})
    add_bridge(bridged_model, LogDetBridge{T})
    add_bridge(bridged_model, RootDetBridge{T})
    add_bridge(bridged_model, RSOCBridge{T})
    add_bridge(bridged_model, RSOCtoPSDBridge{T})
    add_bridge(bridged_model, SOCtoPSDBridge{T})
    return bridged_model
end

include("flip_sign_bridge.jl")
const GreaterToLess{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{GreaterToLessBridge{T}, OT}
const LessToGreater{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{LessToGreaterBridge{T}, OT}
const NonnegToNonpos{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{NonnegToNonposBridge{T}, OT}
const NonposToNonneg{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{NonposToNonnegBridge{T}, OT}
include("vectorizebridge.jl")
const Vectorize{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{VectorizeBridge{T}, OT}
include("scalarizebridge.jl")
const Scalarize{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{ScalarizeBridge{T}, OT}
include("slackbridge.jl")
const ScalarSlack{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{ScalarSlackBridge{T}, OT}
const VectorSlack{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{VectorSlackBridge{T}, OT}
include("functionize_bridge.jl")
const ScalarFunctionize{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{ScalarFunctionizeBridge{T}, OT}
const VectorFunctionize{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{VectorFunctionizeBridge{T}, OT}
include("intervalbridge.jl")
const SplitInterval{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{SplitIntervalBridge{T}, OT}
include("rsocbridge.jl")
const RSOC{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{RSOCBridge{T}, OT}
include("quadtosocbridge.jl")
const QuadtoSOC{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{QuadtoSOCBridge{T}, OT}
include("geomeanbridge.jl")
const GeoMean{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{GeoMeanBridge{T}, OT}
include("square_bridge.jl")
const Square{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{SquareBridge{T}, OT}
include("detbridge.jl")
const LogDet{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{LogDetBridge{T}, OT}
const RootDet{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{RootDetBridge{T}, OT}
include("soctopsdbridge.jl")
const SOCtoPSD{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{SOCtoPSDBridge{T}, OT}
const RSOCtoPSD{T, OT<:MOI.ModelLike} = SingleBridgeOptimizer{RSOCtoPSDBridge{T}, OT}

end # module


text:
# SVG line styles for the different kinds of cut paths of the Shaper Origin:

using DataStructures: OrderedDict

export shaper_cut_attributes, shaper_style_string

shaper_cut_attributes = Dict(
	:inside_cut => OrderedDict(
		"fill" => "white",
		"stroke" => "black",
		"stroke-width" => "0.01",
		"opacity" => "1.0"),
	:outside_cut => OrderedDict(
		"fill" => "black",
		"stroke" => "black",
		"stroke-width" => "0.01",
		"opacity" => "1.0"
	),
	:on_line_cut => OrderedDict(
		"fill" => "none",
		"stroke" => "gray",
		"stroke-width" => "0.01",
		"opacity" => "1.0"
	),
	:on_line_cut => OrderedDict(
		"fill" => "none",
		"stroke" => "gray",
		"stroke-width" => "0.01",
		"opacity" => "1.0"
	),
	:pocket_cut => OrderedDict(
		"fill" => "gray",
		"stroke" => "none",
		"opacity" => "1.0"
	),
	:guide_line => OrderedDict(
		"fill" => "none",
		"stroke" => "blue",
		"stroke-width" => "0.01",
		"opacity" => "1.0"
	)
)

function shaper_style_string(cut_type::Symbol)
    io = IOBuffer()
    for (css, value) in shaper_cut_attributes[cut_type]
        if position(io) > 0
            write(io, " ")
        end
        write(io, "$css: $value;")
    end
    String(take!(io))
end

max_stars_repo_name: MarkNahabedian/ShaperOriginDesignLib

max_issues_repo_name: MarkNahabedian/ShaperOriginDesignLib

max_forks_repo_name: MarkNahabedian/ShaperOriginDesignLib


text:
#=

  Organize a day in Julia ConstraintSolver.jl

  Problem formulation:
  Slides on (finite domain) Constraint Logic Programming, page 38f

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

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
using Printf
const CS = ConstraintSolver
include("constraints_utils.jl")

function organize_day(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    tasks = problem[:tasks]
    durations = problem[:durations]
    precedences = problem[:precedences]
    start_time = problem[:start_time]
    end_time = problem[:end_time]

    n = length(tasks)
    @variable(model, start_time <= begins[1:n] <= end_time, Int)
    @variable(model, start_time <= ends[1:n] <= end_time, Int)

    for i in 1:n
        @constraint(model,ends[i] == begins[i] + durations[i])
    end

    no_overlap(model,begins,durations)

    # precedences
    println("precedences:$precedences")
    for (a,b) in eachrow(precedences)
        @constraint(model, ends[a] <= begins[b])
    end

    @constraint(model,begins[1] >= 11)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                begins_val = convert.(Integer,JuMP.value.(begins; result=sol))
                ends_val = convert.(Integer,JuMP.value.(ends; result=sol))
                # println("begins:$begins_val ")
                # println("ends  :$ends_val ")
                for t in 1:n
                    # println("$(tasks[t]): $(begins_val[t]) .. $(ends_val[t])")
                    @printf "%s: %2d .. %2d\n" tasks[t] begins_val[t] ends_val[t]
                end
                println()

            end
        end
    else
        println("status:$status")
    end

    return status
end

# Problem instance
organize_day_problem = Dict(

    #     task id   1      2      3      4
    :tasks => ["Work","Mail","Shop","Bank"],

    # duration of the four tasks
    :durations => [4,1,2,1],

    # precedences
    # [A,B] : task A must be completed before task B
    :precedences => resize_matrix([[4,3],
                                   [2,1]]
                                  ),
    # Time limits
    :start_time => 9,
    :end_time => 17
)

@time organize_day(organize_day_problem)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# Automatically generated test file
function discontinuous2()

g = [
0.000000000000001,
3.726779962499649,
4.714045207910317,
5.000000000000000,
4.714045207910317,
3.726779962499649,
0.000000000000000,
-5.000000000000000,
-1.000000000000000,
3.000000000000000,
7.000000000000000,
11.000000000000000,
-5.000000000000000,
-1.000000000000000,
3.000000000000000,
7.000000000000000,
11.000000000000000,
0.000000000000001,
3.726779962499649,
4.714045207910317,
5.000000000000000,
4.714045207910317,
3.726779962499649,
0.000000000000000]

ζ_vec = [354.379, 255.565, 184.304, 132.913, 95.8517, 69.1246, 49.85, 35.95, 25.9258, 18.6967]

I_sols = [
[1, 24],
[1, 23, 24],
[1, 14, 16, 24],
[1, 9, 12, 13, 24],
[1, 9, 12, 13, 16, 24],
[1, 5, 8, 12, 13, 16, 24],
[1, 5, 8, 12, 13, 17, 18, 24],
[1, 5, 8, 12, 13, 17, 18, 21, 24]]

f_sols = [
372.169871742722933,
354.378832394845062,
321.422598159772974,
243.754113796361196,
132.333794504652928,
73.691507765404253,
37.521420574787840,
14.150575591039569]

return g, ζ_vec, I_sols, f_sols

end


text:
using Pluto;
Pluto.run();


text:
module JsonGrinder
using Mill, JSON, Printf, Flux
using HierarchicalUtils
import HierarchicalUtils: NodeType, children, InnerNode, LeafNode, printtree, noderepr

using Mill: ArrayNode, BagNode, ProductNode, catobs

const FloatType = Float32

include("schema/schema.jl")
include("extractors/extractors.jl")
include("html_show_tools.jl")
include("hierarchical_utils.jl")
include("util.jl")

export ExtractScalar, ExtractCategorical, ExtractArray, ExtractDict, ExtractOneHot, ExtractVector, MultipleRepresentation, ExtractString, AuxiliaryExtractor, ExtractKeyAsField
export suggestextractor, schema, extractbatch, generate_html

Base.show(io::IO, ::T) where T <: Union{JSONEntry, AbstractExtractor} = show(io, Base.typename(T))
Base.show(io::IO, ::MIME"text/plain", n::Union{JSONEntry, AbstractExtractor}) = HierarchicalUtils.printtree(io, n; trav=false, htrunc=3, vtrunc=20)
Base.getindex(n::Union{JSONEntry, AbstractExtractor}, i::AbstractString) = HierarchicalUtils.walk(n, i)

end # module


text:
using Documenter, GeoPhyInv, SparseArrays, LinearAlgebra, Random, LinearMaps
using Gadfly, Plots, Test

# standard setting for generating doc pages
@init_parallel_stencil(2, false, Float32, 4)

gr()
theme(:juno)

#makedocs(modules=[GeoPhyInv],
#	         doctest=true)
#makedocs()

makedocs(
	 format = Documenter.HTML(
		  prettyurls = get(ENV, "CI", nothing) == "true"),
	 modules=[GeoPhyInv],
   sitename = "GeoPhyInv",
   pages = ["Home" => "index.md",
	    "Medium" => "generated/media/doc.md",
	    "AGeom" => "generated/ageom/doc.md",
	    "SrcWav" => "generated/srcwav/doc.md",
	    "Records" => "generated/data/doc.md",
	    "SeisForwExpt" => Any[
				  "Description" => "generated/fdtd/doc.md",
				  "Tutorial I" => "generated/fdtd/create_snaps.md",
				  "Tutorial II" => "generated/fdtd/reuse_expt.md",
				  ],
	    # "SeisInvExpt" => Any[
		# 		  "Description" => "generated/fwi/doc.md",
		# 		  "Tutorial" => "generated/fwi/pizza.md",

		# 		  "Adjoint Tests" => "generated/fwi/born_tutorial.md",
		# 		 ],
	    # "PoissonExpt" => Any[
				#   "Description" => "generated/Poisson/doc.md",
				#  "Tutorial" => "generated/Poisson/forw.md",
				#  "Adjoint Tests" => "generated/Poisson/test_born.md",
				#  ],
	    ]
#	    "Seismic Born Modeling" => "FWI/born_map.md",
#	    "Seismic Full Waveform Inversion" => "FWI/gradient_accuracy.md",
#    modules = []
)

deploydocs(
    repo   = "github.com/pawbz/GeoPhyInv.jl.git",
)
#=
deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/pawbz/GeoPhyInv.jl.git",
    julia  = "1.0",
    osname = "linux",
    target = "build",
    make   = nothing
)
=#


text:
using VUMPS
using VUMPS:qrpos,lqpos,leftorth,leftenv,rightorth,rightenv,ACenv,Cenv,LRtoC,ALCtoAC,ACCtoALAR,obs_FL,obs_FR,bigleftenv,bigrightenv
using ChainRulesCore
using CUDA
using LinearAlgebra
using OMEinsum
using Random
using Test
using Zygote
CUDA.allowscalar(false)

@testset "matrix autodiff with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64]
    Random.seed!(100)
    a = atype(randn(dtype, 10, 10))
    @test Zygote.gradient(norm, a)[1] ≈ num_grad(norm, a)

    function foo1(x)
        norm(atype(dtype[x 2x; 3x x]))
    end
    @test Zygote.gradient(foo1, 1)[1] ≈ num_grad(foo1, 1)

    # example to solve differential of array of array
    # use `[]` list then reshape
    A = Array{atype,2}(undef, 2, 2)
    for j = 1:2,i = 1:2
        A[i,j] = atype(rand(dtype,2,2))
    end
    function foo2(x)
        # B[i,j] = A[i,j].*x   # mistake
        B = reshape([A[i]*x for i=1:4],2,2)
        return norm(sum(B))
    end
    @test Zygote.gradient(foo2, 1)[1] ≈ num_grad(foo2, 1)
end

@testset "QR factorization with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64]
    Random.seed!(100)
    M = atype(rand(dtype, 3, 3))
    function foo(M)
        Q, R = qrpos(M)
        return norm(Q) + norm(R)
    end
    @test norm(Zygote.gradient(foo, M)[1] - num_grad(foo, M)) ≈ 0 atol = 1e-8
end

@testset "LQ factorization with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64]
    Random.seed!(100)
    M = atype(rand(dtype, 3, 3))
    function foo(M)
        L, Q = lqpos(M)
        return  norm(Q) + norm(L)
    end
    @test norm(Zygote.gradient(foo, M)[1] - num_grad(foo, M)) ≈ 0 atol = 1e-8
end

@testset "loop_einsum mistake with  $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64]
    Random.seed!(100)
    D = 5
    A = atype(rand(dtype, D,D,D))
    B = atype(rand(dtype, D,D))
    function foo(x)
        C = A * x
        D = B * x
        E = ein"abc,abc -> "(C,C)
        F = ein"ab,ab -> "(D,D)
        return norm(Array(E)[]/Array(F)[])
        # E = ein"abc,abc -> "(C,C)[]
        # F = ein"ab,ab -> "(D,D)[]
        # return norm(E/F) mistake for GPU
    end
    @test Zygote.gradient(foo, 1)[1] ≈ num_grad(foo, 1) atol = 1e-8
end

@testset "$(Ni)x$(Nj) leftenv and rightenv with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    D, d = 3, 2
    A = Array{atype{dtype,3},2}(undef, Ni, Nj)
    S = Array{atype{ComplexF64,6},2}(undef, Ni, Nj)
    M = Array{atype{ComplexF64,4},2}(undef, Ni, Nj)
    for j in 1:Nj, i in 1:Ni
        A[i,j] = atype(rand(dtype, D, d, D))
        S[i,j] = atype(rand(ComplexF64, D, d, D, D, d, D))
        M[i,j] = atype(rand(ComplexF64, d, d, d, d))
    end

    ALu, = leftorth(A)
    ALd, = leftorth(A)
    _, ARu, = rightorth(A)
    _, ARd, = rightorth(A)

    function foo1(x)
        _,FL = leftenv(ALu, ALd, M*x)
        s = 0
        for j in 1:Nj, i in 1:Ni
            A = ein"(abc,abcdef),def -> "(FL[i,j], S[i,j], FL[i,j])
            B = ein"abc,abc -> "(FL[i,j], FL[i,j])
            s += norm(Array(A)[]/Array(B)[])
        end
        return s
    end
    @test Zygote.gradient(foo1, 1)[1] ≈ num_grad(foo1, 1) atol = 1e-7

    function foo2(x)
        _,FR = rightenv(ARu, ARd, M*x)
        s = 0
        for j in 1:Nj, i in 1:Ni
            A = ein"(abc,abcdef),def -> "(FR[i,j], S[i,j], FR[i,j])
            B = ein"abc,abc -> "(FR[i,j], FR[i,j])
            s += norm(Array(A)[]/Array(B)[])
        end
        return s
    end
    @test Zygote.gradient(foo2, 1)[1] ≈ num_grad(foo2, 1) atol = 1e-7
end

@testset "$(Ni)x$(Nj) ACenv and Cenv with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    D, d = 3, 2
    A = Array{atype{dtype,3},2}(undef, Ni, Nj)
    S1 = Array{atype{ComplexF64,6},2}(undef, Ni, Nj)
    S2 = Array{atype{ComplexF64,4},2}(undef, Ni, Nj)
    M = Array{atype{ComplexF64,4},2}(undef, Ni, Nj)
    for j in 1:Nj, i in 1:Ni
        A[i,j] = atype(rand(dtype, D, d, D))
        S1[i,j] = atype(rand(ComplexF64, D, d, D, D, d, D))
        S2[i,j] = atype(rand(ComplexF64, D, D, D, D))
        M[i,j] = atype(rand(ComplexF64, d, d, d, d))
    end

    AL, L, _ = leftorth(A)
    R, AR, _ = rightorth(A)
    _, FL = leftenv(AL, AL, M)
    _, FR = rightenv(AR, AR, M)

    C = LRtoC(L, R)
    AC = ALCtoAC(AL, C)
    function foo1(x)
        _, AC = ACenv(AC, FL, M*x, FR)
        s = 0
        for j in 1:Nj, i in 1:Ni
            A = ein"(abc,abcdef),def -> "(AC[i,j], S1[i,j], AC[i,j])
            B = ein"abc,abc -> "(AC[i,j], AC[i,j])
            s += norm(Array(A)[]/Array(B)[])
        end
        return s
    end
    @test Zygote.gradient(foo1, 1)[1] ≈ num_grad(foo1, 1) atol = 1e-8

    function foo2(x)
        _, FL = leftenv(AL, AL, M*x)
        _, FR = rightenv(AR, AR, M*x)
        _, C = Cenv(C, FL, FR)
        s = 0
        for j in 1:Nj, i in 1:Ni
            A = ein"(ab,abcd),cd -> "(C[i,j], S2[i,j], C[i,j])
            B = ein"ab,ab -> "(C[i,j], C[i,j])
            s += norm(Array(A)[]/Array(B)[])
        end
        return s
    end
    @test Zygote.gradient(foo2, 1)[1] ≈ num_grad(foo2, 1) atol = 1e-8
end

@testset "$(Ni)x$(Nj) ACCtoALAR with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    D, d = 3, 2
    A = Array{atype{dtype,3},2}(undef, Ni, Nj)
    S1 = Array{atype{ComplexF64,6},2}(undef, Ni, Nj)
    S2 = Array{atype{ComplexF64,4},2}(undef, Ni, Nj)
    M = Array{atype{ComplexF64,4},2}(undef, Ni, Nj)
    for j in 1:Nj, i in 1:Ni
        A[i,j] = atype(rand(dtype, D, d, D))
        S1[i,j] = atype(rand(ComplexF64, D, d, D, D, d, D))
        S2[i,j] = atype(rand(ComplexF64, D, D, D, D))
        M[i,j] = atype(rand(ComplexF64, d, d, d, d))
    end

    AL, L, _ = leftorth(A)
    R, AR, _ = rightorth(A)
    _, FL = leftenv(AL, AL, M)
    _, FR = rightenv(AR, AR, M)

    Co = LRtoC(L, R)
    ACo = ALCtoAC(AL, Co)
    _, Co = Cenv(Co, FL, FR)
    function foo1(x)
        M *= x
        _, AC = ACenv(ACo, FL, M, FR)
        AL, AR = ACCtoALAR(AC, Co)
        s = 0
        for j in 1:Nj, i in 1:Ni
            A = ein"(abc,abcdef),def -> "(AL[i,j], S1[i,j], AL[i,j])
            B = ein"abc,abc -> "(AL[i,j], AL[i,j])
            s += norm(Array(A)[]/Array(B)[])
            A = ein"(abc,abcdef),def -> "(AR[i,j], S1[i,j], AR[i,j])
            B = ein"abc,abc -> "(AR[i,j], AR[i,j])
            s += norm(Array(A)[]/Array(B)[])
            A = ein"(abc,abcdef),def -> "(AC[i,j], S1[i,j], AC[i,j])
            B = ein"abc,abc -> "(AC[i,j], AC[i,j])
            s += norm(Array(A)[]/Array(B)[])
        end
        return s
    end
    @test Zygote.gradient(foo1, 1)[1] ≈ num_grad(foo1, 1) atol = 1e-2
end

@testset "observable leftenv and rightenv with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    D, d = 3, 2
    A = Array{atype{dtype,3},2}(undef, Ni, Nj)
    S = Array{atype{ComplexF64,6},2}(undef, Ni, Nj)
    M = Array{atype{ComplexF64,4},2}(undef, Ni, Nj)
    for j in 1:Nj, i in 1:Ni
        A[i,j] = atype(rand(dtype, D, d, D))
        S[i,j] = atype(rand(ComplexF64, D, d, D, D, d, D))
        M[i,j] = atype(rand(ComplexF64, d, d, d, d))
    end

    ALu, = leftorth(A)
    ALd, = leftorth(A)
    _, ARu, = rightorth(A)
    _, ARd, = rightorth(A)

    function foo1(x)
        _,FL = obs_FL(ALu, ALd, M*x)
        s = 0
        for j in 1:Nj, i in 1:Ni
            A = ein"(abc,abcdef),def -> "(FL[i,j], S[i,j], FL[i,j])
            B = ein"abc,abc -> "(FL[i,j], FL[i,j])
            s += norm(Array(A)[]/Array(B)[])
        end
        return s
    end
    @test Zygote.gradient(foo1, 1)[1] ≈ num_grad(foo1, 1) atol = 1e-7

    function foo2(x)
        _,FR = obs_FR(ARu, ARd, M*x)
        s = 0
        for j in 1:Nj, i in 1:Ni
            A = ein"(abc,abcdef),def -> "(FR[i,j], S[i,j], FR[i,j])
            B = ein"abc,abc -> "(FR[i,j], FR[i,j])
            s += norm(Array(A)[]/Array(B)[])
        end
        return s
    end
    @test Zygote.gradient(foo2, 1)[1] ≈ num_grad(foo2, 1) atol = 1e-7
end

@testset "bigleftenv and bigrightenv with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    D, d = 3, 2
    A = Array{atype{dtype,3},2}(undef, Ni, Nj)
    S = Array{atype{ComplexF64,8},2}(undef, Ni, Nj)
    M = Array{atype{ComplexF64,4},2}(undef, Ni, Nj)
    for j in 1:Nj, i in 1:Ni
        A[i,j] = atype(rand(dtype, D, d, D))
        S[i,j] = atype(rand(ComplexF64, D, d, d, D, D, d, d, D))
        M[i,j] = atype(rand(ComplexF64, d, d, d, d))
    end

    ALu, = leftorth(A)
    ALd, = leftorth(A)
    _, ARu, = rightorth(A)
    _, ARd, = rightorth(A)
    function foo1(x)
        _, FL4 = bigleftenv(ALu, ALd, M*x)
        s = 0
        for j in 1:Nj, i in 1:Ni
            A = ein"(abcd,abcdefgh),efgh -> "(FL4[i,j], S[i,j], FL4[i,j])
            B = ein"abcd,abcd -> "(FL4[i,j], FL4[i,j])
            s += norm(Array(A)[]/Array(B)[])
        end
        return s
    end
    @test Zygote.gradient(foo1, 1)[1] ≈ num_grad(foo1, 1) atol = 1e-7

    function foo2(x)
        _, FR4 = bigrightenv(ARu, ARd, M*x)
        s = 0
        for j in 1:Nj, i in 1:Ni
            A = ein"(abcd,abcdefgh),efgh -> "(FR4[i,j], S[i,j], FR4[i,j])
            B = ein"abcd,abcd -> "(FR4[i,j], FR4[i,j])
            s += norm(Array(A)[]/Array(B)[])
        end
        return s
    end
    @test Zygote.gradient(foo2, 1)[1] ≈ num_grad(foo2, 1) atol = 1e-7
end


text:
module test_locales

include("../../src/locales.jl")
using .locales
using Test
locale_map = locales._locale_map

function test_locale_data_structure()
        for locale_cls in values(locales._locale_map)
            @test length(locale_cls.language_tag) == 2
            @test all(islowercase, locale_cls.language_tag)
            @test length(locale_cls.month_names) == 12
            @test all(!isempty, locale_cls.month_names)  # not blank or None
#            @test all(!is(nothing), locale_cls.month_names)  # not blank or None
            @test length(locale_cls.gregorian_month_names) == 12
 #           @test all(locale_cls.gregorian_month_names)  # not blank or None
            @test length(locale_cls.day_names) == 7
  #          @test all(locale_cls.day_names)  # not blank or None
            @test locale_cls.notation != nothing
            @test locale_cls.gregorian_notation != nothing
	end
end

function test_locale_map()
        @test length(locales._locale_map) > 0
        @test "en" in keys(locales._locale_map)
        @test "ar" in keys(locales._locale_map)
        @test "az" in keys(locales._locale_map)
        @test "bn" in keys(locales._locale_map)
end
end

#function test_duplicated_language_tag()
 #       @test_throws LookupError
#
 #           class ExtraLocale(locales.Locale):
  #              language_tag = "en"

max_stars_repo_name: ceferisbarov/hijri-converter

max_issues_repo_name: ceferisbarov/hijri-converter


text:
# plots behavioral data

using Printf, Gadfly, Colors
import Cairo, Fontconfig

include("behaviorstats.jl")

const nrats = 1   # treat all 4 rats as one rat
const t = Theme(minor_label_font_size=6pt, major_label_font_size=8pt)
const figpath = "../fitdata/fits"

function plotdataset(dataset::Symbol, stimx, magctrst, stimlabel, ctrstlabel, ttitle, ratid::Int = 0)
    d = load_data(dataset, ratid)
    # compute stats
    d_bystim, d_bymag = getpsychchronstats(d, nrats)
    d_seq, d_seqfit, d_biasid, d_biasmag = getseqstats(d, nrats)
    # generate plots
    pchronid, ppsychid = plotchronpsychid(d_bystim, stimx, stimlabel, ttitle, t)
    pchronmag, ppsychmag = plotchronpsychmag(d_bymag, magctrst, ctrstlabel, ttitle, t)
    pcondpsych = plotcondpsych(d_seq, d_seqfit, stimx, stimlabel, ttitle, t)
    pbiasid = plotbiasid(d_biasid, stimx, "prev. $stimlabel", ttitle, t)
    pbiasmag = plotbiasmag(d_biasmag, magctrst, "prev. $ctrstlabel", ttitle, t)

    # plot all of those onto single "sheet"
    p = vstack(
        hstack(vstack(pchronid, ppsychid), vstack(pchronmag, ppsychmag)),
        vstack(
            hstack(pcondpsych[1], pcondpsych[2]),
            hstack(pcondpsych[3], pcondpsych[4])
        ),
        hstack(pbiasid, pbiasmag)
    )
    # output to file
    outfile = ratid == 0 ? @sprintf("%s/behav_%s.pdf", figpath, string(dataset)) :
                           @sprintf("%s/behav_%s_rat%d.pdf", figpath, string(dataset), ratid)
    draw(PDF(outfile, 8.5inch, 11inch), p)
end

function plotitlvddataset(ratid::Int = 0)
    mixctrst = [0.12, 0.36, 0.6, 1.0]
    nconcentr = 4
    d = load_data(:itlvd, ratid)
    # generate plot layers, iterating over concetrations
    lchron = Array{Gadfly.Layer}(undef, nconcentr, 3)
    lpsych = Array{Gadfly.Layer}(undef, nconcentr, 3)
    for i = 1:nconcentr  # iterate over different concentrations
        # stats for given concentration
        _, d_bymag = getpsychchronstats(d[d[:concentr] .== i,:], nrats)
        d_bymag[:tmin] = d_bymag[:tμ] - d_bymag[:tsem]
        d_bymag[:tmax] = d_bymag[:tμ] + d_bymag[:tsem]
        d_bymag[:pcorrmin] = d_bymag[:pcorrμ] - d_bymag[:pcorrsem]
        d_bymag[:pcorrmax] = d_bymag[:pcorrμ] + d_bymag[:pcorrsem]
        # the 8i-36 turns the stimmag into the (1-4) range
        d_bymag[:stimctrst] = mixctrst[d_bymag[:stimmag] .+ (8i - 36)]
        ti = copy(t)
        ti.default_color = RGB(([1.0, 1.0, 1.0] .- 0.2i)...)
        lchron[i,1] = layer(d_bymag, x="stimctrst", y="tμ", Geom.point, ti)[1]
        lchron[i,2] = layer(d_bymag, x="stimctrst", y="tμ", ymin="tmin", ymax="tmax",
                            Geom.errorbar, ti)[1]
        lchron[i,3] = layer(d_bymag, x="stimctrst", y="tμ", Geom.line, ti)[1]
        lpsych[i,1] = layer(d_bymag, x="stimctrst", y="pcorrμ", Geom.point, ti)[1]
        lpsych[i,2] = layer(d_bymag, x="stimctrst", y="pcorrμ", ymin="pcorrmin", ymax="pcorrmax",
                            Geom.errorbar, ti)[1]
        lpsych[i,3] = layer(d_bymag, x="stimctrst", y="pcorrμ", Geom.line, ti)[1]
    end
    # assemble plots
    pchronid = plot(lchron[:,1], lchron[:,2], lchron[:,3],
        Guide.xlabel("concentration contrast", orientation=:horizontal),
        Guide.ylabel("reaction time [s]", orientation=:vertical),
        Guide.title("Chronometric curve, interleaved"),
        Scale.x_continuous(minvalue=0.0, maxvalue=1.0),
        Scale.y_continuous(minvalue=0.26, maxvalue=0.44))
    ppsychid = plot(lpsych[:,1], lpsych[:,2], lpsych[:,3],
        Guide.xlabel("concentration contrast", orientation=:horizontal),
        Guide.ylabel("correct choiced (fraction)", orientation=:vertical),
        Guide.title("Psychometric curve, interleaved"),
        Scale.x_continuous(minvalue=0.0, maxvalue=1.0),
        Scale.y_continuous(minvalue=0.5, maxvalue=1.0))
    p = vstack(pchronid, ppsychid)
    draw(PDF("$(figpath)/behav_itlvd.pdf", 8.5inch, 11inch), p)
end

const categ_x = [0, 0.2, 0.32, 0.44, 0.56, 0.68, 0.8, 1]
#const categ_ctrst = categ_x - flipdim(categ_x, 1)
const categ_ctrst = [0, 0, 0, 0, 0.12, 0.36, 0.6, 1.0]
plotdataset(:categ, categ_x, categ_ctrst,
            "odor A (frac)", "concentration contrast", "categorization")
for ratid in [1, 2, 3, 4]
    plotdataset(:categ, categ_x, categ_ctrst,
                "odor A (frac)", "concentration contrast",
                "categorization, rat$(ratid)", ratid)
end

const ident_x = [-4, -3, -2, -1, 1, 2, 3, 4]
const ident_ctrst = [-5.0, -5, -5, -5, -4, -3, -2, -1]
plotdataset(:ident, ident_x, ident_ctrst,
            "concentration (10^x)", "concentration contrast", "identification")
for ratid in [1, 2, 3, 4]
    plotdataset(:ident, ident_x, ident_ctrst,
                "concentration (10^x)", "concentration contrast",
                "identification, rat$(ratid)", ratid)
end

plotitlvddataset()


text:
"""
Euler functions for entropy variables and more

"""

# 1D wavespeed for use in interface fluxes
function wavespeed(rho,rhou,E)
    cvel = (@. sqrt(γ*pfun(rho,rhou,E)/rho))
    return (@. sqrt(abs(rhou/rho)) + cvel)
end

Unorm(U) = sum((x->x.^2).(U))

"primitive pressure to conservative vars"
function primitive_to_conservative(rho,U,p)
    rhoU = (x->rho.*x).(U)
    unorm = Unorm(U)
    E = (@. p/(γ-1) + .5*rho*unorm)
    return (rho,rhoU...,E)
end

function primitive_to_conservative(rho,u,v,p)
    return primitive_to_conservative(rho,(u,v),p)
end
function primitive_to_conservative(rho,u,v,w,p)
    return primitive_to_conservative(rho,(u,v,w),p)
end

"inverse temperature (used in entropy conservative fluxes)"
function betafun(rho,rhoU,E)
    p = pfun(rho,rhoU,E)
    return (@. rho/(2*p))
end
function betafun(rho,rhou,rhov,E)
    return betafun(rho,(rhou,rhov),E)
end
function betafun(rho,rhou,rhov,rhow,E)
    return betafun(rho,(rhou,rhov,rhow),E)
end

"pressure as a function of ρ,u,v,E"
function pfun(rho,rhoU,E,rhounorm)
    return @. (γ-1)*(E-.5*rhounorm)
end
function pfun(rho,rhoU,E)
    rhounorm = Unorm(rhoU)./rho
    return pfun(rho,rhoU,E,rhounorm)
end

# hack
function pressure_fun(rho,rhou,rhov,E)
    rhounorm = @. (rhou^2+rhov^2)/rho
    return @. (γ-1)*(E-.5*rhounorm)
    #return pfun(rho,(rhou,rhov),E)
end

"specific energy as a function of conservative variables"
function rhoefun(rho,rhoU,E)
    rhoUnorm = Unorm(rhoU)
    return (@. E - .5*rhoUnorm/rho)
end

"Thermodynamic entropy as a function of conservative variables"
function sfun(rho,rhoU,E)
    rhoe = rhoefun(rho,rhoU,E)
    return (@. log((γ-1)*rhoe/(rho^γ)))
end

"Mathematical entropy"
function Sfun(rho,rhou,rhov,E)
    return -rho.*sfun(rho,(rhou,rhov),E)
end
function Sfun(rho,rhou,rhov,rhow,E)
    return -rho.*sfun(rho,(rhou,rhov,rhow),E)
end

"Entropy variables as functions of conservative vars"
function v_ufun(rho,rhoU,E)
    ρe = rhoefun(rho,rhoU,E)
    sU = sfun(rho,rhoU,E)
    v1 = (@. (-E + ρe*(γ + 1 - sU))/ρe)
    vU = [@. rhoUi/ρe for rhoUi in rhoU]
    vE = (@. (-rho)/ρe)
    return (v1,vU...,vE)
end
function v_ufun(rho,rhou,rhov,E)
    return v_ufun(rho,(rhou,rhov),E)
end
function v_ufun(rho,rhou,rhov,rhow,E)
    return v_ufun(rho,(rhou,rhov,rhow),E)
end

"entropy as function of entropy variables"
function s_vfun(v1,vU,vE)
    vUnorm = Unorm(vU)
    return (@. γ - v1 + vUnorm/(2*vE))
end

"specific energy as function of entropy variables"
function rhoe_vfun(v1,vU,vE)
    s = s_vfun(v1,vU,vE)
    return (@. ((γ-1)/((-vE)^γ))^(1/(γ-1)) * exp(-s/(γ-1)))
end

"Conservative vars as functions of entropy variables"
function u_vfun(v1,vU,vE)
    rhoeV = rhoe_vfun(v1,vU,vE)
    vUnorm = Unorm(vU)
    rho   = (@. rhoeV*(-vE))
    rhoU  = [@. rhoeV*(vUi) for vUi in vU]
    E     = (@. rhoeV*(1-vUnorm/(2*vE)))
    return (rho,rhoU...,E)
end
function u_vfun(v1,v2,v3,v4)
    return u_vfun(v1,(v2,v3),v4)
end
function u_vfun(v1,v2,v3,v4,v5)
    return u_vfun(v1,(v2,v3,v4),v5)
end

max_stars_repo_name: jlchan/ESDG-CNS

max_issues_repo_name: jlchan/ESDG-CNS

max_forks_repo_name: jlchan/ESDG-CNS


text:
### mcvar and mcse stand for Monte Carlo variance and Monte Carlo error respectively

### Monte Carlo variance (erroneously) assuming IID samples

mcvar_iid(x::Vector{Float64}) = var(x)/length(x)

mcvar_iid(x::Matrix{Float64}, pars::Range=1:size(x, 2)) = Float64[mcvar_iid(x[:, pars[i]]) for i = 1:length(pars)]

mcvar_iid(x::Matrix{Float64}, par::Real) = mcvar_iid(x, par:par)

mcvar_iid(c::MCChain, pars::Range=1:size(c.samples, 2)) = mvcvar_iid(c.samples, pars)

mcvar_iid(c::MCChain, par::Real) = mcvar_iid(c, par:par)

### Monte Carlo standard error (erroneously) assuming IID samples

msce_iid(x::Vector{Float64}) = sqrt(mcvar_iid(x))

msce_iid(x::Matrix{Float64}, pars::Range=1:size(x, 2)) = Float64[mcse_iid(x[:, pars[i]]) for i = 1:length(pars)]

msce_iid(x::Matrix{Float64}, par::Real) = msce_iid(x, par:par)

msce_iid(c::MCChain, pars::Range=1:size(c.samples, 2)) = msce_iid(c.samples, pars)

msce_iid(c::MCChain, par::Real) = msce_iid(c, par:par)

### Function for estimating Monte Carlo variance using batch means, see for instance Flegal J.M, Jones, G.L. Batch Means
### and Spectral Variance Estimators in Markov chain Monte Carlo. Annals of Statistics, 2010, 38 (2), pp 1034-1070

function mcvar_bm(x::Vector{Float64}; batchlen::Int=100)
  nbatches = div(length(x), batchlen)
  @assert nbatches > 1 "Choose batch size such that the number of batches is greather than one"
  nbsamples = nbatches*batchlen
  batchmeans = Float64[mean(x[((j-1)*batchlen+1):(j*batchlen)]) for j = 1:nbatches]
  return batchlen*var(batchmeans)/nbsamples
end

mcvar_bm(x::Matrix{Float64}, pars::Range=1:size(x, 2); batchlen::Int=100) =
  Float64[mcvar_bm(x[:, pars[i]]; batchlen=batchlen) for i = 1:length(pars)]

mcvar_bm(x::Matrix{Float64}, par::Real; batchlen::Int=100) = mcvar_bm(x, par:par; batchlen=batchlen)

mcvar_bm(c::MCChain, pars::Range=1:size(c.samples, 2); batchlen::Int=100) =
  mcvar_bm(c.samples, pars; batchlen=batchlen)

mcvar_bm(c::MCChain, par::Real; batchlen::Int=100) = mcvar_bm(c, par:par; batchlen=batchlen)

### Monte Carlo standard error using batch means

mcse_bm(x::Vector{Float64}; batchlen::Int=100) = sqrt(mcvar_bm(x; batchlen=batchlen))

mcse_bm(x::Matrix{Float64}, pars::Range=1:size(x, 2); batchlen::Int=100) =
  Float64[mcse_bm(x[:, pars[i]]; batchlen=batchlen) for i = 1:length(pars)]

mcse_bm(x::Matrix{Float64}, par::Real; batchlen::Int=100) = mcse_bm(x, par:par; batchlen=batchlen)

mcse_bm(c::MCChain, pars::Range=1:size(c.samples, 2); batchlen::Int=100) =
  mcse_bm(c.samples, pars; batchlen=batchlen)

mcse_bm(c::MCChain, par::Real; batchlen::Int=100) = mcse_bm(c, par:par; batchlen=batchlen)

### Function for estimating the variance of a single MC chain using the initial monotone sequence estimator (IMSE), see
### Geyer C.J. Practical Markov Chain Monte Carlo. Statistical Science, 1992, 7 (4), pp 473-483

function mcvar_imse(x::Vector{Float64}; maxlag::Int=length(x)-1)
  k = convert(Int, floor((maxlag-1)/2))
  m = k+1

  # Preallocate memory
  g = Array(Float64, k+1)

  # Calculate empirical autocovariance
  acv = autocov(x, 0:maxlag)

  # Calculate \hat{G}_{n, m} from the autocovariances, see pp. 477 in Geyer
  for j = 0:k
    g[j+1] = acv[2*j+1]+acv[2*j+2]
    if g[j+1] <= 0
      m = j
      break
    end
  end

  # Create the monotone sequence of g
  if m > 1
    for j = 2:m
      if g[j] > g[j-1]
        g[j] = g[j-1]
      end
    end
  end

  # Calculate the initial monotone sequence estimator
  return (-acv[1]+2*sum(g[1:m]))/length(x)
end

mcvar_imse(x::Matrix{Float64}, pars::Range=1:size(x, 2); maxlag::Int=size(x, 1)-1) =
  Float64[mcvar_imse(x[:, pars[i]]; maxlag=maxlag) for i = 1:length(pars)]

mcvar_imse(x::Matrix{Float64}, par::Real; maxlag::Int=size(x, 1)-1) = mcvar_imse(x, par:par; maxlag=maxlag)

mcvar_imse(c::MCChain, pars::Range=1:size(c.samples, 2); maxlag::Int=size(c.samples, 1)-1) =
  mcvar_imse(c.samples, pars; maxlag=maxlag)

mcvar_imse(c::MCChain, par::Real; maxlag::Int=size(c.samples, 1)-1) = mcvar_imse(c, par:par; maxlag=maxlag)

### Monte Carlo standard error using Geyer's IMSE

mcse_imse(x::Vector{Float64}; maxlag::Int=length(x)-1) = sqrt(mcvar_imse(x; maxlag=maxlag))

mcse_imse(x::Matrix{Float64}, pars::Range=1:size(x, 2); maxlag::Int=size(x, 1)-1) =
  Float64[mcse_imse(x[:, pars[i]]; maxlag=maxlag) for i = 1:length(pars)]

mcse_imse(c::MCChain, pars::Range=1:size(c.samples, 2); maxlag::Int=size(c.samples, 1)-1) =
  mcse_imse(c.samples, pars; maxlag=maxlag)

mcse_imse(c::MCChain, par::Real; maxlag::Int=size(c.samples, 1)-1) = mcse_imse(c, par:par; maxlag=maxlag)

### Function for estimating the variance of a single MC chain using the initial positive sequence estimator (IPSE), see
### Geyer C.J. Practical Markov Chain Monte Carlo. Statistical Science, 1992, 7 (4), pp 473-483

function mcvar_ipse(x::Vector{Float64}; maxlag::Int=length(x)-1)
  k = convert(Int, floor((maxlag-1)/2))

  # Preallocate memory
  g = Array(Float64, k+1)
  m = k+1

  # Calculate empirical autocovariance
  acv = autocov(x, 0:maxlag)

  # Calculate \hat{G}_{n, m} from the autocovariances, see pp. 477 in Geyer
  for j = 0:k
    g[j+1] = acv[2*j+1]+acv[2*j+2]
    if g[j+1] <= 0
      m = j
      break
    end
  end

  # Calculate the initial monotone sequence estimator
  return (-acv[1]+2*sum(g[1:m]))/length(x)
end

mcvar_ipse(x::Matrix{Float64}, pars::Range=1:size(x, 2); maxlag::Int=size(x, 1)-1) =
  Float64[mcvar_ipse(x[:, pars[i]]; maxlag=maxlag) for i = 1:length(pars)]

mcvar_ipse(x::Matrix{Float64}, par::Real; maxlag::Int=size(x, 1)-1) = mcvar_ipse(x, par:par; maxlag=maxlag)

mcvar_ipse(c::MCChain, pars::Range=1:size(c.samples, 2); maxlag::Int=size(c.samples, 1)-1) =
  mcvar_ipse(c.samples, pars; maxlag=maxlag)

mcvar_ipse(c::MCChain, par::Real; maxlag::Int=size(c.samples, 1)-1) = mcvar_ipse(c, par:par; maxlag=maxlag)

### Monte Carlo standard error using Geyer's IMSE

mcse_ipse(x::Vector{Float64}; maxlag::Int=length(x)-1) = sqrt(mcvar_ipse(x; maxlag=maxlag))

mcse_ipse(x::Matrix{Float64}, pars::Range=1:size(x, 2); maxlag::Int=size(x, 1)-1) =
  Float64[mcse_ipse(x[:, pars[i]]; maxlag=maxlag) for i = 1:length(pars)]

mcse_ipse(c::MCChain, pars::Range=1:size(c.samples, 2); maxlag::Int=size(c.samples, 1)-1) =
  mcse_ipse(c.samples, pars; maxlag=maxlag)

mcse_ipse(c::MCChain, par::Real; maxlag::Int=size(c.samples, 1)-1) = mcse_ipse(c, par:par; maxlag=maxlag)

### Wrapper function for computing MC variance using various approaches

vtypes = (:bm, :iid, :imse, :ipse)

function mcvar(x::Vector{Float64}; vtype::Symbol=:imse, args...)
  @assert in(vtype, vtypes) "Unknown variance type $vtype"

  if vtype == :bm
    return mcvar_bm(x; args...)
  elseif vtype == :iid
    return mcvar_iid(x)
  elseif vtype == :imse
    return mcvar_imse(x; args...)
  elseif vtype == :ipse
    return mcvar_ipse(x; args...)
  end
end

mcvar(x::Matrix{Float64}, pars::Range=1:size(x, 2); vtype::Symbol=:imse, args...) =
  Float64[mcvar(x[:, pars[i]]; vtype=vtype, args...) for i = 1:length(pars)]

mcvar(x::Matrix{Float64}, par::Real; vtype::Symbol=:imse, args...) = mcvar(x, par:par; vtype=vtype, args...)

mcvar(c::MCChain, pars::Range=1:size(c.samples, 2); vtype::Symbol=:imse, args...) =
  mcvar(c.samples, pars; vtype=vtype, args...)

mcvar(c::MCChain, par::Real; vtype::Symbol=:imse, args...) = mcvar(c, par:par; vtype=vtype, args...)

### Wrapper function for computing Monte Carlo (standard) error using various approaches

function mcse(x::Vector{Float64}; vtype::Symbol=:imse, args...)
  @assert in(vtype, vtypes) "Unknown stdiance type $vtype"

  if vtype == :bm
    return mcse_bm(x; args...)
  elseif vtype == :iid
    return mcse_iid(x)
  elseif vtype == :imse
    return mcse_imse(x; args...)
  elseif vtype == :ipse
    return mcse_ipse(x; args...)
  end
end

mcse(x::Matrix{Float64}, pars::Range=1:size(x, 2); vtype::Symbol=:imse, args...) =
  Float64[mcse(x[:, pars[i]]; vtype=vtype, args...) for i = 1:length(pars)]

mcse(x::Matrix{Float64}, par::Real; vtype::Symbol=:imse, args...) = mcse(x, par:par; vtype=vtype, args...)

mcse(c::MCChain, pars::Range=1:size(c.samples, 2); vtype::Symbol=:imse, args...) =
  mcse(c.samples, pars; vtype=vtype, args...)

mcse(c::MCChain, par::Real; vtype::Symbol=:imse, args...) = mcse(c, par:par; vtype=vtype, args...)


text:
# Read configration files

###################################
"""
    params = geodata("simlation/path/_output"::AbstractString)
    params = geodata("simlation/path/_output/geoclaw.data"::AbstractString)

Function: geoclaw.data reader
"""
function geodata(outdir::AbstractString)
    ## set filename
    fname = occursin("geoclaw.data", basename(outdir)) ? outdir : joinpath(outdir, "geoclaw.data")

    ## check whether it exists
    if !isfile(fname); error("File $fname is not found."); end
    ## read all lines
    open(fname,"r") do f
        global txt = readlines(f)
    end
    ## parse parameters
    # parameters (mandatory?)
    cs = parse(Float64,split(txt[occursin.("coordinate",txt)][1],r"\s+")[1])
    p0 = parse(Float64,split(txt[occursin.("ambient_pressure",txt)][1],r"\s+")[1])
    R = parse(Float64,split(txt[occursin.("earth_radius",txt)][1],r"\s+")[1])
    eta0 = parse(Float64,split(txt[occursin.("sea_level",txt)][1],r"\s+")[1])
    dmin = parse(Float64,split(txt[occursin.("dry_tolerance",txt)][1],r"\s+")[1])
    # parameters (optional?)
    if any(occursin.("manning_coefficient",txt))
        n = parse(Float64,split(txt[occursin.("manning_coefficient",txt)][1],r"\s+")[1])
    else
        n = 0.0
    end
    ## instance
    params = VisClaw.GeoParam(cs,p0,R,eta0,n,dmin)
    ## return values
    return params
end
###################################

###################################
"""
    surgeparams = surgedata("simlation/path/_output"::AbstractString)
    surgeparams = surgedata("simlation/path/_output/surge.data"::AbstractString)

Function: surge.data reader
"""
function surgedata(outdir::AbstractString)
    ## set filename
    fname = occursin("surge.data", basename(outdir)) ? outdir : joinpath(outdir, "surge.data")

    ## check whether it exists
    if !isfile(fname); error("File $fname is not found."); end
    ## read all lines
    open(fname,"r") do f
        global txt = readlines(f)
    end
    ## parse parameters
    windindex = parse(Int64,split(txt[occursin.("wind_index",txt)][1],r"\s+")[1])
    slpindex = parse(Int64,split(txt[occursin.("pressure_index",txt)][1],r"\s+")[1])
    stormtype = parse(Int64,split(txt[occursin.("storm_specification_type",txt)][1],r"\s+")[1])
    ## instance
    surgeparams = VisClaw.SurgeParam(windindex,slpindex,stormtype)
    ## return values
    return surgeparams
end
###################################

###################################
"""
    gaugeinfo = gaugedata("simlation/path/_output"::AbstractString)
    gaugeinfo = gaugedata("simlation/path/_output/gauges.data"::AbstractString)

Function: gauges.data reader
"""
function gaugedata(outdir::AbstractString)
    ## set filename
    fname = occursin("gauges.data", basename(outdir)) ? outdir : joinpath(outdir,"gauges.data")

    ## check whether it exists
    if !isfile(fname); error("File $fname is not found."); end
    ## read all lines
    open(fname,"r") do f
        global txt = readlines(f)
    end
    ## parse parameters
    ngauges = parse(Int64, split(txt[occursin.("ngauges",txt)][1],r"\s+")[1])
    ## preallocate
    gaugeinfo = Vector{VisClaw.Gauge}(undef,ngauges)

    ## read gauge info
    baseline = findfirst(x->occursin("ngauges", x), txt)
    for i = 1:ngauges
        txtline = split(strip(txt[baseline+i]), r"\s+", keepempty=false)
        label = txtline[1]
        id = parse(Int64,txtline[1])
        loc = [parse(Float64,txtline[2]), parse(Float64,txtline[3])]
        time = [parse(Float64,txtline[4]), parse(Float64,txtline[5])]
        # instance
        gaugeinfo[i] = VisClaw.Gauge(label,id,0,loc,[],time,[])
    end

    ## return values
    return gaugeinfo
end
###################################

###################################
"""
    amrparam = amrdata("simlation/path/_output"::AbstractString)
    amrparam = amrdata("simlation/path/_output/amr.data"::AbstractString)

Function: amr.data reader
"""
function amrdata(outdir::AbstractString)
    ## set filename
    fname = occursin("amr.data", basename(outdir)) ? outdir : joinpath(outdir, "amr.data")

    ## check whether it exists
    if !isfile(fname); error("File $fname is not found."); end
    ## read all lines
    open(fname,"r") do f
        global txt = readlines(f)
    end

    ## parse parameters
    maxlevel = parse(Int64,split(txt[occursin.("amr_levels_max",txt)][1],r"\s+")[1])
    ## instance
    amrparam = VisClaw.AMRParam(maxlevel)
    ## return values
    return amrparam
end
###################################

###################################
"""
    regions = regiondata("simlation/path/_output"::AbstractString)
    regions = regiondata("simlation/path/_output/region.data"::AbstractString)

Function: region.data reader
"""
function regiondata(outdir::AbstractString)
    ## set filename
    fname = occursin("regions.data", basename(outdir)) ? outdir : joinpath(outdir, "regions.data")

    ## check whether it exists
    if !isfile(fname); error("File $fname is not found."); end
    ## read all lines
    open(fname,"r") do f
        global txt = readlines(f)
    end

    ## parse parameters
    nregion = parse(Int64, split(txt[occursin.("num_regions",txt)][1],r"\s+")[1])
    ## preallocate
    regions = Vector{VisClaw.Region}(undef,nregion)

    ## read gauge info
    baseline = findfirst(x->occursin("num_regions", x), txt)
    for i = 1:nregion
        txtline = split(strip(txt[baseline+i]), r"\s+", keepempty=false)
        minlevel = parse(Int64,txtline[1])
        maxlevel = parse(Int64,txtline[2])
        tl = (parse(Float64,txtline[3]), parse(Float64,txtline[4]))
        xl = (parse(Float64,txtline[5]), parse(Float64,txtline[6]))
        yl = (parse(Float64,txtline[7]), parse(Float64,txtline[8]))
        ## instance
        regions[i] = VisClaw.Region(minlevel:maxlevel, tl, xl, yl)
    end

    return regions
end
###################################


text:
using RobustAndOptimalControl, ControlSystems

Gtrue = tf([11.2], [1, 0.12, 0])
M, wB, A = 1.5, 20.0, 1e-8
WS = tf([1 / M, wB], [1, wB * A])

WU = ss(1.0)

WT = []

P_test = let
    _A = [-2.0e-7 0.0 -22.4; 0.0 -0.12 0.0; 0.0 0.25 0.0]
    _B = [4.0 0.0; 0.0 8.0; 0.0 0.0]
    _C = [4.99999996666667 0.0 -3.73333333333333; 0.0 0.0 0.0; 0.0 0.0 -5.6]
    _D = [0.666666666666667 0.0; 0.0 1.0; 1.0 0.0]
    ss(_A, _B, _C, _D)
end

P = RobustAndOptimalControl.hinfpartition(Gtrue, WS, WU, WT)
@test tf(ss(P)) ≈ tf(P_test)

P = RobustAndOptimalControl.hinfpartition2(Gtrue, WS, WU, WT)
@test tf(ss(P)) ≈ tf(P_test)

## Flexible servo example

A = [
    -0.015948101119319734 0.0 0.04185849126302005 0.0
    0.0 -0.011069870050970518 0.0 0.03334113387614947
    0.0 0.0 -0.04185849126302005 0.0
    0.0 0.0 0.0 -0.03334113387614947
]
B = [0.08325 0.0; 0.0 0.0628125; 0.0 0.04785714285714286; 0.031218750000000007 0.0]
C = [0.5 0.0 0.0 0.0; 0.0 0.5 0.0 0.0]
D = [0.0 0.0; 0.0 0.0]
G = ss(A, B, C, D)

A = [-0.001 0.0; 0.0 -0.001]
B = [1.0 0.0; 0.0 1.0]
C = [0.0999 0.0; 0.0 0.0999]
D = [0.1 0.0; 0.0 0.1]
WS = ss(A, B, C, D)

D = [0.01 0.0; 0.0 0.01]
WU = ss(D)

A = [-1.0 0.0; 0.0 -1.0]
B = [1.0 0.0; 0.0 1.0]
C = [-9.9 0.0; 0.0 -9.9]
D = [10.0 0.0; 0.0 10.0]
WT = ss(A, B, C, D)

P_test = let
    _A = [
        -0.001 0.0 0.0 0.0 -0.5 0.0 0.0 0.0
        0.0 -0.001 0.0 0.0 0.0 -0.5 0.0 0.0
        0.0 0.0 -1.0 0.0 0.5 0.0 0.0 0.0
        0.0 0.0 0.0 -1.0 0.0 0.5 0.0 0.0
        0.0 0.0 0.0 0.0 -0.0159481011193197 0.0 0.04185849126302 0.0
        0.0 0.0 0.0 0.0 0.0 -0.0110698700509705 0.0 0.0333411338761495
        0.0 0.0 0.0 0.0 0.0 0.0 -0.04185849126302 0.0
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0333411338761495
    ]
    _B = [
        1.0 0.0 0.0 0.0
        0.0 1.0 0.0 0.0
        0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0
        0.0 0.0 0.08325 0.0
        0.0 0.0 0.0 0.0628125
        0.0 0.0 0.0 0.0478571428571429
        0.0 0.0 0.03121875 0.0
    ]
    _C = [
        0.0999 0.0 0.0 0.0 -0.05 0.0 0.0 0.0
        0.0 0.0999 0.0 0.0 0.0 -0.05 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 -9.9 0.0 5.0 0.0 0.0 0.0
        0.0 0.0 0.0 -9.9 0.0 5.0 0.0 0.0
        0.0 0.0 0.0 0.0 -0.5 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 -0.5 0.0 0.0
    ]
    _D = [
        0.1 0.0 0.0 0.0
        0.0 0.1 0.0 0.0
        0.0 0.0 0.01 0.0
        0.0 0.0 0.0 0.01
        0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0
        1.0 0.0 0.0 0.0
        0.0 1.0 0.0 0.0
    ]
    ss(_A, _B, _C, _D)
end

P = RobustAndOptimalControl.hinfpartition(G, WS, WU, WT)
@test_broken tf(ss(P)) ≈ tf(P_test)

P = RobustAndOptimalControl.hinfpartition2(G, WS, WU, WT)
@test tf(ss(P)) ≈ tf(P_test)


text:
using Brownian
using Base.Test

p = BrownianMotion([0, 0.1, 0.5, 0.75, 1])

rand(p)
rand([p, p])


text:
# to create vector, put comma
V = [1, 2, 3, 4, 5, 6]

# to create matrix just add space
V1 = [1 2 3 4 5 6]

# to declare an array just put semicolon
M1 = [1 2 3; 4 5 6; 7 8 9]

# to resahepe vector or matrix just use resahepe()

V3 = [1, 2, 3, 4, 5, 6, 7, 8, 9]
reshape(V3, 3,3) # for vector

M2 = [1 2 3; 4 5 6]
reshape(M2, 3, 2) # for matrix 2x3 to 3x2

# to create identitiy matirx
using LinearAlgebra
Matrix{Int}(I, 4, 4)

# for dimensions and size of matrix use ndims() and size(), and to see type of elements in matrix using eltype()
ndims(M2)
size(M2)
eltype(M1)

# to concatenate two vector use vcat (similar to cbind() in R) and use hcat similar to (rbind() in R)

d1 = [1 2 3; 4 5 6]
d2 = [6 7 8; 9 10 11]

hcat(d1, d2)
vcat(d1, d2)

# functions in julia

function f(x, y, z)
    x + y + z
end

# or just assign

g(a, b, c) = a + b + c

# function with return

function β(a, b, c)
    result = a + b + c
    return result
end

β(2,π,3)

# we can use opetator as function

2 * 3 * 4
*(2,3,4)

# function with default parameters

z(x, y, a = 1, b = 2) = a*x + b*y
z(3,5)

 # Anonymous function

x -> x + 1 # or

function (x)
    x +1
end

# Anonymous function generally useing with other first class function
map(x -> x + 1, [1, 2])

# multiple return values

function ff( a, b, c)
    a * b * c, a + b + c
end;

ff(1, 2, 3)

x, y  = ff(1, 2, 3);

display(x)
display(y)

function fff(a, b, c)
    multi = a * b * c
    add = a + b + c

    return multi, add
end

fff(1, 2, 3)

# using Anonymous function with multiple piece of code, do-blocks will help

map(x -> x + 1, [1, 2, 3, 4]) # basic way

map([1, 2, 3, 4]) do x # other way
     x + 1
end

# control flow
# compound expression

x = (a = 5; b = 6; a*b)

# other way

x = begin
    a = 5
    b = 6
    a * b
end

# conditional evaluation

function ineq(x, y)
    if x > y
        x = 2
    elseif x < y
        x = 1
    else
        x = 0
    end
end

ineq(1, 2)

# repeated evaluatins: loops

j = 3;

while j > 0
    println(j^2)
    global  j -= 1
end

# for loops

for j in [1, 2, 3]
    println(j)
end

# we can add ∈ to loops

for j ∈ [1, 2, 3]
    println(j)
end

# warning and informational messages

function f(x)
    if x<1
        @warn "x must be positive!"
        else
            return x^2
    end
end

f(2)

function g(x)
    if x < 1
        @info "X must be positive!"
    else
        return x^2
    end
end

g(-2)

# Plotting

using PyPlot
x = 1:100
y = rand(100)
display(gcf())
plot(x,y, color = "tomato")
title("basic graphs in PyPlot")

max_stars_repo_name: fmyilmaz/DSGE

max_issues_repo_name: fmyilmaz/DSGE

max_forks_repo_name: fmyilmaz/DSGE


text:
import PWS
using PWS: Trajectory, SMCEstimate, DirectMCEstimate, marginal_configuration, MarginalEnsemble, ConditionalEnsemble, gene_expression_system, generate_configuration, logpdf, cooperative_chemotaxis_system, reaction_network