text:
#=

Draw radar map in resulting world frame (after solve)

=#

##

using Distributed
# addprocs(10)

using Colors
using Caesar
@everywhere using Caesar

using DistributedFactorGraphs: listDataEntrySequence

import Rotations as _Rot
using CoordinateTransformations

# using DistributedFactorGraphs.DocStringExtensions
using Dates
using JSON2
# using BSON
using Serialization
using FixedPointNumbers
using Manifolds
using StaticArrays
using TensorCast
using DataStructures

using ImageMagick, FileIO
using Images
using ImageDraw
using ImageView
using ImageTransformations
using ImageFeatures

using VideoIO

using RoMEPlotting
Gadfly.set_default_plot_size(35cm,20cm)

using Optim

## load the solved graph and reattach the datastores

dfg_datafolder = "/tmp/caesar/philos"

fg = loadDFG("$dfg_datafolder/results_5/x0_5_x280")

##

ds = FolderStore{Vector{UInt8}}(:radar, "$dfg_datafolder/data/radar")
addBlobStore!(fg, ds)

ds = FolderStore{Vector{UInt8}}(:gps_fix, "$dfg_datafolder/data/gps")
addBlobStore!(fg, ds)

ds = FolderStore{Vector{UInt8}}(:lidar, "$dfg_datafolder/data/lidar")
addBlobStore!(fg, ds)

ds = FolderStore{Vector{UInt8}}(:camera, "$dfg_datafolder/data/camera")
addBlobStore!(fg, ds)

## structure

# walk through all the poses in fg
#  pull out all world frame radar images
#  pull out all camera images at each (not all poses have world location solution)
#  draw side by side radar on left and camera on right (more camera than radar images)
#

##

M = SpecialEuclidean(2)
e0 = identity_element(M)

##

# get two radar sweep point clouds and rotate them to the world frame
function _fetchDataPointCloudWorld(fg::AbstractDFG, lb::Symbol)
  de,db = getData(fg, lb, :RADAR_SWEEP)
  # deserialize is not legacy safe, rather use BSON.@save(io, Caesar._PCL.PCLPointCloud2(pc))
  bPC = deserialize(PipeBuffer(db))
  wPb = exp(M, e0, hat(M, e0, getPPE(fg, lb).suggested))
  apply(M, wPb, bPC)
end

Base.clamp(img::AbstractMatrix{<:RGB},b::Real,t::Real) = map(px->RGB{N0f8}(clamp(px.r,b,t),clamp(px.g,b,t),clamp(px.b,b,t)), img)

function buildPoseImages( fg::AbstractDFG,
                          lb::Symbol;
                          # _X_ = 3000,
                          # _Y_ = 1200,
                          pattern::Regex=r"IMG_CENTER_",
                          T = RGB{N0f8},
                          color=Gray{Float64}(0.05),
                          scale_map=1.2,
                          ppose=zeros(3) )
  #
  global img_wPC, imgL

  M = SpecialEuclidean(2)
  e0 = identity_element(M)

  R_ = _Rot.RotMatrix(0.0) # 0.017*pi

  framestack = []

  # fetch the camera images
  ents = listDataEntrySequence(fg, lb, pattern, sortDFG)
  imgs_lb = fetchDataImage.(fg, lb, ents)
  rows,cols = size(imgs_lb[1])

  # expand the radar map
  if isSolved(fg, lb)
    wPC__ = _fetchDataPointCloudWorld(fg, lb);
    # add traj plotting here (draw points and lines on radar plot)

    nRw = _Rot.RotMatrix(140*pi/180.0)
    nTw = ProductRepr([0.0;0.0], nRw)
    nPC = apply(M, nTw, wPC__)
    __ppe = getPPE(fg, lb).suggested
    _ppe = Manifolds.compose(M, nTw, exp(M, e0, hat(M, e0, __ppe)))
    ppe = vee(M, e0, log(M, e0, _ppe))
    img_wPC_ = makeImage!(nPC, (-2500,1000),(-2000,1000); rows=trunc(Int, scale_map*rows), color, pose=ppe, ppose=ppose );
    ppose .= ppe
    if lb == :x0
      img_wPC = deepcopy(img_wPC_)
    else
      img_wPC .+= img_wPC_
    end
    imgG = (px->RGB(0,10*px,0)).(img_wPC_);

    imgL = RGB.(img_wPC) + imgG
  end

  imgs_R = map(img->warp(img, ImageTransformations.recenter(R_, center(img))), imgs_lb)

  frame_ = Matrix{T}(undef, size(imgs_R[1])...)

  # write mosaic frames to the video
  for frame in imgs_R
    ax, ay = axes(frame)
    I = ax.parent .+ ax.offset
    J = ay.parent .+ ay.offset
    for (i,i_) in enumerate(I), (j,j_) in enumerate(J)
      frame_[i,j] = frame[i_,j_]
    end
    imgF = mosaic([imgL, frame_]; nrow=1) #, npad=50)
    push!(framestack, imgF)
  end

  # typeof(framestack[1])
  framestack_ = Vector{Matrix{T}}(undef, length(framestack));
  for (i,fr) in enumerate(framestack)
    framestack_[i] = clamp(fr, 0,1);
  end

  return framestack_
end

##

lbls = sortDFG(ls(fg)) #[1:41]

# imgL = nothing    # Matrix{RGB{Float64}}()
global img_wPC = nothing
global imgL = nothing

ppose = zeros(3)
pattern = r"STAB_IMG_CENTER_"

framestack = buildPoseImages(fg, :x0; ppose, pattern );

##

encoder_options = (crf=23, preset="medium")
framerate=90
open_video_out("/tmp/caesar/philos/mosaic_stab.mp4", framestack[1], framerate=framerate, encoder_options=encoder_options) do writer

  for (i,lb) in enumerate(lbls)
    println("doing all frames for ", lb)

    framestack = buildPoseImages(fg, lb; ppose, pattern )

    for frame in framestack
      # frame_ .= T.(frame)
      write(writer, frame)
    end

  end # for

end # video

##

lb = :x48
des = listDataEntrySequence(fg, lb, r"IMG_CENTER", sortDFG)
imgs = fetchDataImage.(fg, lb, des)

##

#


text:
#!/home/jburdick/bin/julia
# Writes out the coefficients of the bound matrix, both
# using exact binomials, and using various approximations
# (and log-transformed numbers), as a check of the math.

push!(LOAD_PATH, ".")

using CountingBound
using LogCountingBound

"""
  Writes counts for some values of k, r, and n.
  k: the value of k
  maxN: the maximum value of n
  outputDir: directory in which to write output files
  Side effects: writes files A and b.
    (The choices of r are a bit arbitrary).
"""
function writeCounts(k, maxN, outputDir)
  # values of n to use, from k to maxN
  # (for now, spaced at powers of 2)
  nList = [BigInt(2^i) for i in 1:trunc(log(2,maxN))]
  print(nList)
  print("\n")
  print(typeof(nList[1]))
  print("\n")
  # ??? for some reason, this was giving a syntax error
  # nList = [n for n in nList if n>=k && n<=maxN]
  # therefore, using filter()...
  nList = filter(n -> n>=k && n<=maxN, nList)

  # FIXME should create output directory
  # first, write the bound b (and its log)
  of = open(outputDir * "/b_k=" * string(k) * "_maxN=" * string(maxN) * ".csv", "w")
  write(of, "k,n,bound,exp(log(bound)),log(bound)\n")
  for n in nList
    bound = countingBound(binomial(n, 2), binomial(n, k))
		logBound = NaN
		if bound > 0
			logBound = log(bound)
		end
    write(of, string(k) * "," * string(n) * "," * string(bound) * ","
			* string(exp(logBound)) * ","
      * string(logBound) * "\n")
  end
  close(of)

  # then, write the coefficients
  of = open(outputDir * "/A_k=" * string(k) * "_maxN=" * string(maxN) * ".csv", "w")
  write(of, "k,r,n,A,exp(logA),logA\n")
  for n = nList
    # FIXME what should the bound on r be?
    for r = k:min(n, 2*k)
      print("k=" * string(k) * " r=" * string(r) * " n=" * string(n) * "\n")
      A = approxNumMaximalCliques1(k, r, n)
      logA = logNumMaximalCliques(k, r, n)
      write(of, string(k) * "," * string(r) * "," * string(n)
        * "," * string(A)
				* "," * string(exp(logA))
        * "," * string(logA) * "\n")
    end
  end
  close(of)

end

mkpath("coef1")
# smoke tests
writeCounts(2, 200, "coef1")
writeCounts(3, 200, "coef1")
# almost a smoke test
# writeCounts(4, 100, "coef1")
writeCounts(8, 1000, "coef1")
writeCounts(10, 2000, "coef1")
writeCounts(12, 4000, "coef1")

ext: jl

max_stars_repo_name: joshtburdick/misc

max_issues_repo_name: joshtburdick/misc

max_forks_repo_name: joshtburdick/misc


text:
#!/usr/bin/env julia

using GZip

# Function for reading input .vcf file
function read(file)
	if ismatch(r".vcf$", file)   	# Normal .vcf files
		f = open(file)
		lines = readlines(f)
		close(f)
	elseif ismatch(r".gz$", file)   # Gzipped .vcf files
		f = GZip.open(file)
		lines = readlines(f)
		close(f)
	end
	return lines
end


text:
n_points = 10
tau = 1

@testset "Refining triangulation" begin
    embedding = embedding_example(n_points, 3, tau)

    triang_points, simplex_indices = triangulate(embedding)

    # How many simplices are there in the triangulation?
    n_orig = size(simplex_indices, 1)
    simplices = [triang_points[simplex_indices[i, :], :] for i in 1:n_orig]

    @testset "Triangulation" begin
        # Create an 3D embedding to triangulate.
        @test isa(triang_points, Array{Float64, 2})
        @test isa(simplex_indices, Array{Int, 2})
    end

    # Rules for splitting simplices in 3D with size reducing factor 2
    splitting_rules = simplicial_subdivision(2, 3)

    @testset "Simplex division rules" begin
        @test isa(splitting_rules[1], Array{Int64, 2})
        @test isa(splitting_rules[2], Array{Int64, 2})
    end

    @testset "Refine triangulation" begin
        @testset "E = $E" for E in 2:4
            @testset "k = $k" for k in 1:3
                splitting_rules = simplicial_subdivision(k, E)
                embedding = embedding_example(n_points, E, 1)
                triang_points, simplex_indices = triangulate(embedding)

                # How many simplices did we get?
                n_orig = size(simplex_indices, 1)

                simplices = [triang_points[simplex_indices[i, :], :] for i in 1:n_orig]
                refined_triangulation = refine_triangulation(triang_points, simplex_indices, collect(1:n_orig), k)

                # How many new simplices did we get?
                n_refined = size(refined_triangulation[2], 1)
                refined_simplices = [refined_triangulation[1][refined_triangulation[2][i, :], :]
                                        for i in 1:n_refined]

                refined_vertices = refined_triangulation[1]
                refined_simplex_indices = refined_triangulation[2]

                @testset "Arrays are of right dimensions and types" begin
                    @test isa(refined_vertices, Array{Float64, 2})
                    @test isa(refined_simplex_indices, Array{Int, 2})
                end
                # Compare volume of original triangulation with the refined triangulation

                @testset "Volume is preserved" begin
                    triang_vols = [abs(det(hcat(simplices[i], ones(E + 1)))) for i in 1:n_orig]
                    refined_triang_vols = [abs(det(hcat(refined_simplices[i], ones(E + 1)))) for i in 1:n_refined]
                    @test sum(triang_vols) ≈ sum(refined_triang_vols) atol = 1/10^7
                end
            end
        end
    end

    @testset "Simplex volumes are split by factor of k" begin

        @testset "k = $k" for k in 1:3
            @testset "E = $E" for E in 2:7
                canonical_simplex_vertices = zeros(E + 1, E)
                canonical_simplex_vertices[2:(E + 1), :] = eye(E)

                # Refine the canonical simplex
                simplex_indices = zeros(Int, 1, E + 1)
                simplex_indices[1, :] = collect(1:E+1)

                # Generic splitting rules for this combination of dimension and splitting
                # factor.
                splitting_rules = simplicial_subdivision(k, E)

                refined = refine_triangulation(canonical_simplex_vertices, simplex_indices, [1], k)

                # How many new simplices did we get?
                n = size(refined[2], 1)
                refined_simplices = [refined[1][refined[2][i, :], :] for i in 1:n]

                canonical_simplex_vol = abs(det(
                    hcat(canonical_simplex_vertices, ones(E + 1))
                    ))

                volumes_subsimplices = [abs(det(hcat(refined_simplices[i], ones(E + 1))))
                                                                            for i in 1:n]

                @test canonical_simplex_vol ≈ sum(volumes_subsimplices) atol = 1/10^7
            end
        end
    end
end

@testset "Centroids and radii of triangulation" begin
    @testset "E = $E" for E in 2:4
        @testset "k = $k" for k in 1:4
            splitting_rules = simplicial_subdivision(k, E)
            embedding = embedding_example(n_points, E, 1)
            triang_points, simplex_indices = triangulate(embedding)

            centroids, radii = centroids_radii2(triang_points, simplex_indices)
            @test isa(centroids, Array{Float64, 2})
            @test isa(radii, Array{Float64, 1})
        end
    end
end

max_stars_repo_name: kahaaga/SimplexSplitting

max_issues_repo_name: kahaaga/SimplexSplitting

max_forks_repo_name: kahaaga/SimplexSplitting


text:
# ===============================
# Written by AAE
# TU Delft, Spring 2014
# simulkade.com
# ===============================

# ================================================================
# Changes:
#    2015-01-10 changed numberofcells to dims
# ================================================================

# ======================= Linear source term ========================
function linearSourceTerm(betta0::CellValue)
m = betta0.domain
d = m.dimension
if (d ==1) || (d==1.5)
  Nx = m.dims[1]
  G = [1:Nx+2;]
  b = betta0.value[2:end-1]
  row_index = reshape(G[2:Nx+1],Nx)  # main diagonal (only internal cells)
  AP_diag = reshape(b,Nx)
  sparse(row_index, row_index, AP_diag, Nx+2, Nx+2)
elseif (d == 2) || (d == 2.5) || (d==2.8)
  Nx = m.dims[1]
  Ny = m.dims[2]
  G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)
  b = betta0.value[2:end-1,2:end-1]
  row_index = reshape(G[2:Nx+1,2:Ny+1],Nx*Ny)  # main diagonal (only internal cells)
  AP_diag = reshape(b,Nx*Ny)
  sparse(row_index, row_index, AP_diag, (Nx+2)*(Ny+2), (Nx+2)*(Ny+2))
elseif (d == 3) || (d==3.2)
  Nx = m.dims[1]
  Ny = m.dims[2]
  Nz = m.dims[3]
  G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)
  b = betta0.value[2:end-1,2:end-1,2:end-1]
  row_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],Nx*Ny*Nz)  # main diagonal (only internal cells)
  AP_diag = reshape(b,Nx*Ny*Nz)
  sparse(row_index, row_index, AP_diag, (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
end
end

function linearSourceTerm(m::MeshStructure, betta0::Real)
linearSourceTerm(CellValue(m, betta0))
end

function linearSourceTerm{T<:Real}(m::MeshStructure, betta0::Array{T})
d = m.dimension
if (d ==1) || (d==1.5)
  Nx = m.dims
  G = [1:Nx+2;]
  row_index = reshape(G[2:Nx+1],Nx)  # main diagonal (only internal cells)
  AP_diag = reshape(betta0,Nx)
  sparse(row_index, row_index, AP_diag, Nx+2, Nx+2)
elseif (d == 2) || (d == 2.5) || (d==2.8)
  Nx = m.dims[1]
  Ny = m.dims[2]
  G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)
  row_index = reshape(G[2:Nx+1,2:Ny+1],Nx*Ny)  # main diagonal (only internal cells)
  AP_diag = reshape(betta0,Nx*Ny)
  sparse(row_index, row_index, AP_diag, (Nx+2)*(Ny+2), (Nx+2)*(Ny+2))
elseif (d == 3) || (d==3.2)
  Nx = m.dims[1]
  Ny = m.dims[2]
  Nz = m.dims[3]
  G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)
  row_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],Nx*Ny*Nz)  # main diagonal (only internal cells)
  AP_diag = reshape(betta0,Nx*Ny*Nz)
  sparse(row_index, row_index, AP_diag, (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
end
end

# ================================== constant source term ================================
function constantSourceTerm(phi0::CellValue)
m = phi0.domain
d = m.dimension
if (d ==1) || (d==1.5)
  Nx = m.dims[1]
  G = [1:Nx+2;]
  row_index = reshape(G[2:Nx+1],Nx)  # main diagonal (only internal cells)
  RHS = zeros(Nx+2)
  RHS[row_index] = reshape(phi0.value[2:end-1],Nx)
elseif (d == 2) || (d == 2.5) || (d==2.8)
  Nx = m.dims[1]
  Ny = m.dims[2]
  G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)
  row_index = reshape(G[2:Nx+1,2:Ny+1],Nx*Ny)  # main diagonal (only internal cells)
  RHS = zeros((Nx+2)*(Ny+2))
  RHS[row_index] = reshape(phi0.value[2:end-1,2:end-1],Nx*Ny)
elseif (d == 3) || (d==3.2)
  Nx = m.dims[1]
  Ny = m.dims[2]
  Nz = m.dims[3]
  G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)
  row_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],Nx*Ny*Nz)  # main diagonal (only internal cells)
  RHS = zeros((Nx+2)*(Ny+2)*(Nz+2))
  RHS[row_index] = reshape(phi0.value[2:end-1,2:end-1,2:end-1],Nx*Ny*Nz)
end
RHS
end

function constantSourceTerm(m::MeshStructure, phi0::Real)
constantSourceTerm(CellValue(m, phi0))
end

function constantSourceTerm{T<:Real}(m::MeshStructure, phi0::Array{T})
d = m.dimension
if (d ==1) || (d==1.5)
  Nx = m.dims
  G = [1:Nx+2;]
  row_index = reshape(G[2:Nx+1],Nx)  # main diagonal (only internal cells)
  RHS = zeros(Nx+2)
  RHS[row_index] = reshape(phi0,Nx)
elseif (d == 2) || (d == 2.5) || (d==2.8)
  Nx = m.dims[1]
  Ny = m.dims[2]
  G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)
  row_index = reshape(G[2:Nx+1,2:Ny+1],Nx*Ny)  # main diagonal (only internal cells)
  RHS = zeros((Nx+2)*(Ny+2))
  RHS[row_index] = reshape(phi0,Nx*Ny)
elseif (d == 3) || (d==3.2)
  Nx = m.dims[1]
  Ny = m.dims[2]
  Nz = m.dims[3]
  G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)
  row_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],Nx*Ny*Nz)  # main diagonal (only internal cells)
  RHS = zeros((Nx+2)*(Ny+2)*(Nz+2))
  RHS[row_index] = reshape(phi0,Nx*Ny*Nz)
end
RHS
end


text:
### saving parameters

function saveparameters(path_to_results, T, number_zones, cars_per_zone, C, e_drive, p_min, p_max, e_dest, A_drive)

    # preparing the parameter values for saving
    parameters = zeros(1, 9);

    parameters[1,1] = T;
    parameters[1,2] = number_zones;
    parameters[1,3] = cars_per_zone;
    parameters[1,4] = C;
    parameters[1,5] = e_drive;
    parameters[1,6] = p_min;
    parameters[1,7] = p_max;
    parameters[1,8] = e_dest;
    parameters[1,9] = A_drive;

    parameters = convert(DataFrame, parameters);

    # create the header for the parameter file and save the results
    header_vector = ["T (time steps)", "number_zones", "cars_per_zone", "C (number of cars)", "e_drive (model parameter 1)", "p_min (model parameter 2)", "p_max (model parameter 3)", "e_dest (model parameter 4)", "A_drive" ];
    results_string = string(path_to_results,"/sampling_parameters.csv");
    CSV.write(results_string, parameters, header=header_vector);

end

max_stars_repo_name: ArsamAryandoust/CarParkingMaps

max_issues_repo_name: ArsamAryandoust/CarParkingMaps

max_forks_repo_name: ArsamAryandoust/CarParkingMaps


text:
function test_importexport(::Type{T}) where T <: CasadiSymbolicObject
    @testset "$( string("Import export functions ", T, "                       ") )" begin
        M = rand(3,3)
        R = rand(5,3)

        @test to_julia( inv( T(M) ) )  ≈ inv(M)
        @test to_julia( sqrt( T(M) ) ) ≈ sqrt.(M)
        @test to_julia( sin( T(M) ) )  ≈ sin.(M)
        @test to_julia( cos( T(M) ) )  ≈ cos.(M)
        @test to_julia( vec( T(M) ) )  ≈ vec(M)
        @test to_julia(T(R)')  ≈ R'

        @test size( T(M) ) == size(M)
        @test to_julia( reshape( T(M), (9,1) ) ) ≈ reshape( M, (9,1) )
    end
end


text:
using SpreadsheetFunctions
using DataFrames

sheet3_table = DataFrame([[1, 2, 3, 4, 5], [6.7, 8.2, 3.4, 5.6, 2.1], ["a", "b", "c", "d", "e"]], [:a, :b, :c])
sheet3_table.d_amt = map(row -> row.a+row.b, eachrow(sheet3_table))
sheet3_e2 = SUM(sheet3_table[2, [:a, :b, :c, :d_amt]])
sheet3_e4 = PRODUCT(sheet3_table[2:3, [:a, :b]])
sheet3_a7 = SUM(sheet3_table.a)
sheet3_b7 = SUM(sheet3_table.b)
sheet3_c7 = SUM(sheet3_table.c)
sheet3_d7 = SUM(sheet3_table.d_amt)


text:
using Random
using Functors

include("field_grid.jl")
include("radfuncs.jl")
Random.seed!(1)

mutable struct LinearOperator
    ranks::Any
    rmin::AbstractFloat
    rmax::AbstractFloat

    radfunc::Any
    conv::Function

    filter::AbstractArray
    grid::Grid
end
@functor LinearOperator

"""
    function LinearOperator(
        name;
        dx = nothing,
        rmax = nothing,
        ranks = nothing,
        grid = nothing,
        dims = 3,
        radfunc = nothing,
        rmin = 0.0,
        σ = 1.0
    )
"""
function LinearOperator(
    name;
    dx = nothing,
    rmax = nothing,
    ranks = nothing,
    grid = nothing,
    dims = 3,
    radfunc = nothing,
    rmin = 0.0,
    σ = 1.0,
)
name=Symbol(name)
    if name == :neural
        radfunc = Radfunc(; rmin, rmax)
    elseif name == :potential
        li = lo = lf = 0
        rmin = 2dx
        radfunc = r-> 1 / (4π * r)
    elseif name == :inverse_squared_field
        li = 0
        lf = lo = 1
        rmin = 2dx
        radfunc = r-> 1 / (4π * r^2)
    elseif name == :Gaussian
        li = lf = lo = 1
        radfunc = r-> exp(-r^2 / (2σ^2)) / (2π)^(dims / 2)
    elseif name == :grad
        dV = dx^dims
        li = 0
        lf = lo = 1
        rmax = dx
        n = 1
        radfunc = r-> δn(r, dx, n, dims)
    elseif name == :div
        dV = dx^dims
        lf = li = 1
        lo = 0
        rmax = dx
        n = 1
        radfunc = r-> δn(r, dx, n, dims)
    elseif name == :curl
        dV = dx^dims
        lf = lo = li = 1
        rmax = dx
        n = 1
        radfunc = r-> -δn(r, dx, n, dims)
    elseif name == :Laplacian
        dV = dx^dims
        lf = li = lo = 0
        rmax = dx
        n = 2
        radfunc = r-> δn(r, dx, n, dims)
    end

    if ranks === nothing
        ranks = (li, lf, lo)
    end
    if grid === nothing
        grid = Grid(dx, rmax; dims, rank_max = ranks[2])
    end
    conv = FieldConv(ranks...)

if typeof(radfunc)===Radfunc
    radfunc1=radfunc
else
    radfunc1(r) = (rmin - tol) < r < (rmax + tol) ? radfunc(r) : 0.0
end

    filter = Field(; radfunc=radfunc1, grid, rank = ranks[2])
     LinearOperator(ranks, rmin, rmax, radfunc1, conv, filter, grid)
end

"""
    function (m::LinearOperator)(x::AbstractArray, grid::Grid)

"""
function (m::LinearOperator)(x::AbstractArray)
    if m.rmax == 0
        return m.radfunc(0) * x
    end

    ix = 0
    x_ = 0
    dV = 0.0
    Zygote.ignore() do
        ix = [
            Int(a):Int(b + a - 1) for (a, b) in zip(m.grid.origin, size(x))
            # Int(a):Int(b + a - 1) for (a, b) in zip(y.grid.origin, x.grid.sz)
        ]
        x_ = x
        # x_ = x.components
        dV = m.grid.dV
    end
    y_ = m.filter

    res = m.conv(x_, y_)
    # res* dV
    res[ix...,:] .* dV
end
# makefilter(radfunc,grid) = Field(;radfunc, grid)
function remake!(m)
    rank, radfunc, grid = field_rank(m.filter), m.radfunc, m.grid
    m.filter = Field(; rank, radfunc, grid)
end

max_stars_repo_name: aced-differentiate/equiv-docs

max_issues_repo_name: aced-differentiate/equiv-docs

max_forks_repo_name: aced-differentiate/equiv-docs


text:
#=
Copyright 2020 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

"""
    runSimulation(input_filename::String; verbose::Bool=false, out_files::Bool=false)

Execute the simulation main loop.

Args:
    - `input_filename`: The name of the `.yml` input file.
    - `verbose`: Opt. Boolean flag for STDOUT. Default is `false`.
    - `out_files`: Opt. Boolean flag to control the `.out` files writing. Default `false`.
"""
function runSimulation(input_filename::String; verbose::Bool=false, out_files::Bool=false)
    data = loadSimulationFiles(input_filename)
    blood = buildBlood(data["blood"])

    verbose && println("Build $input_filename arterial network \n")

    jump = data["solver"]["jump"]

    vessels, edges = buildArterialNetwork(data["network"], blood, jump)
    makeResultsFolder(data, input_filename)

    Ccfl = data["solver"]["Ccfl"]
    heart = vessels[1].heart
    total_time = data["solver"]["cycles"]*heart.cardiac_T
    timepoints = range(0, stop=heart.cardiac_T, length=jump)

    verbose && println("Start simulation \n")

    current_time = 0.0
    passed_cycles = 0

    verbose && (@printf("Solving cardiac cycle no: %02d", passed_cycles + 1); starting_time = time_ns())

    counter = 1
    conv_criteria = "norm"
    conv_toll = data["solver"]["convergence tolerance"]

    while true
        dt = calculateDeltaT(vessels, Ccfl)
        solveModel(vessels, edges, blood, dt, current_time)
        updateGhostCells(vessels)

        if current_time >= timepoints[counter]
            saveTempData(current_time, vessels, counter)
            counter += 1
        end

        if (current_time - heart.cardiac_T*passed_cycles) >= heart.cardiac_T &&
          (current_time - heart.cardiac_T*passed_cycles + dt) > heart.cardiac_T

            if passed_cycles + 1 > 1
                err, err_loc = computeConvError(vessels)
                verbose && printConvError(err, err_loc)
            else
                print('\n')
            end

            transferTempToLast(vessels)

            out_files && transferLastToOut(vessels)

            if (passed_cycles+1>1 && checkConvergence(err, conv_toll))
	        writeConv(data, passed_cycles)
		break
            end

            passed_cycles += 1
            verbose && @printf("Solving cardiac cycle no: %02d", passed_cycles + 1)

            timepoints = timepoints .+ heart.cardiac_T
            counter = 1
        end

        current_time += dt
        if current_time >= total_time
            passed_cycles += 1
            verbose && println("\nNot converged after $passed_cycles cycles, End!")
            break
        end
    end
    verbose && (@printf("\n"); ending_time = (time_ns() - starting_time)/1.0e9)
    verbose && println("Elapsed time = $ending_time seconds")

    writeResults(vessels)

    cd("..")
end


text:
function pid(::Type{T}, stimulus::AbstractVector{Int}, responses::AbstractMatrix{Int}) where T
    pid!(Hasse(T, size(responses,1)), stimulus, responses; zero=false)
end

function pid(::Type{T}, stimulus::AbstractVector{Int},
             responses::AbstractMatrix{Int},
             names::AbstractVector) where T
    if length(names) != size(responses,1)
        throw(ArgumentError("number of names provided does not match the number of responses"))
    end

    pid!(Hasse(T, names), stimulus, responses; zero=false)
end

function pid!(h::Hasse, df::DataFrame, stimulus, responses, args...; kwargs...)
    pid!(h, df[:, stimulus], transpose(Array(df[:, responses])), args...; kwargs...)
end

function pid(::Type{T}, df::DataFrame, stimulus, responses, args...; kwargs...) where T
    pid(T, df[:, stimulus], transpose(Array(df[:, responses])), responses, args...; kwargs...)
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown]
# https://github.com/JuliaStats/StatsModels.jl/issues/220

# %% [markdown]
# ## Reproduce the slow-down (1)

# %%
using DataFrames
using GLM
using StatsBase: sample

n = 20
r = 10
x_symbols = [Symbol("x$i") for i in 1:n]
@time df = DataFrame(rand(100, n+1), [:y; x_symbols]);

# %%
x_vars = sample(x_symbols, r; replace=false)
@time F = term(:y) ~ sum(term(x) for x in x_vars)
@time cols = Tables.columntable(df)
@time mf = ModelFrame(F, cols, model=LinearModel)
@time mm = ModelMatrix(mf)
@time y = response(mf)
@time linmodel = fit(LinearModel, mm.m, y)
@time regmodel = StatsModels.TableRegressionModel(linmodel, mf, mm)

# %%
x_vars = sample(x_symbols, r; replace=false)
@time F = term(:y) ~ sum(term(x) for x in x_vars)
@time cols = Tables.columntable(df)
@time mf = ModelFrame(F, cols, model=LinearModel)
@time mm = ModelMatrix(mf)
@time y = response(mf)
@time linmodel = fit(LinearModel, mm.m, y)
@time regmodel = StatsModels.TableRegressionModel(linmodel, mf, mm)

# %% [markdown]
# ## Reproduce the slow-down (2)

# %%
using DataFrames
using GLM
using StatsBase: sample

n = 20
r = 10
x_symbols = [Symbol("x$i") for i in 1:n]
df = DataFrame(rand(100, n+1), [:y; x_symbols]);

result = []
for _ in 1:10
    x_vars = sample(x_symbols, r; replace=false)
    F = term(:y) ~ sum(term(x) for x in x_vars)
    @time regmodel = lm(F, df)
    push!(result, regmodel)
end
result

# %% [markdown]
# ## Solution

# %%
using DataFrames
using GLM
using StatsBase: sample

struct MyLinearModel{T, Y, X}
    linmodel::T
    y_var::Y
    x_vars::X
end

function my_lm(y_var::Symbol, x_vars::AbstractVector{Symbol}, df::DataFrame)
    y = df[!, y_var]
    X = [ones(nrow(df)) Matrix(df[!, x_vars])]
    linmodel = lm(X, y)
    MyLinearModel(linmodel, y_var, x_vars)
end

get_y_var(F::FormulaTerm) = F.lhs.sym
get_x_vars(F::FormulaTerm) = collect((t -> t.sym).(F.rhs))
my_lm(F::FormulaTerm, df::DataFrame) = my_lm(get_y_var(F), get_x_vars(F), df)

function Base.show(io::IO, mylm::MyLinearModel)
    linmodel, y_var, x_vars = getfield.(Ref(mylm), (:linmodel, :y_var, :x_vars))
    ct = coeftable(linmodel)
    ct.rownms .= string.((Symbol("(Intercept)"), x_vars...,))
    print(io, typeof(mylm), "\n\n")
    print(io, y_var, " ~ 1")
    for x in x_vars print(io, " + ", x) end
    print(io, "\n\nCoefficients:\n")
    show(io, ct)
    print(io, "\n")
end

n = 20
r = 10
x_symbols = [Symbol("x$i") for i in 1:n]
df = DataFrame(rand(100, n+1), [:y; x_symbols]);

myresult = []
for _ in 1:10
    x_vars = sample(x_symbols, r; replace=false)
    F = term(:y) ~ sum(term(x) for x in x_vars)
    @time mylinmodel = my_lm(F, df)
    push!(myresult, mylinmodel)
end
myresult

# %% [markdown]
# ## Comparison

# %%
F = @formula(y ~ x9 + x5 + x1 + x4 + x6 + x8 + x7 + x2 + x3)
@time lm(F, df)

# %%
F = @formula(y ~ x9 + x5 + x1 + x4 + x6 + x8 + x7 + x3 + x2)
@time my_lm(F, df)

# %%
F = @formula(y ~ x9 + x5 + x1 + x4 + x6 + x8 + x3 + x7 + x2)
@time my_lm(F, df)

# %%
F = @formula(y ~ x9 + x5 + x1 + x4 + x6 + x8 + x3 + x2 + x7)
y_var, x_vars = get_y_var(F), get_x_vars(F)
@time my_lm(y_var, x_vars, df)

# %% [markdown]
# ## Analysis of the slow-down

# %%
x_vars = sample(x_symbols, r; replace=false)
@show x_vars
@time F = term(:y) ~ sum(term(x) for x in x_vars)
@time lm(F, df)

# %%
@time my_lm(F, df)

# %%
x_vars = sample(x_symbols, r; replace=false)
@time F = term(:y) ~ sum(term(x) for x in x_vars)
@time my_lm(F, df)

# %%
x_vars = sample(x_symbols, r; replace=false)
@time F = term(:y) ~ sum(term(x) for x in x_vars)
@time cols = Tables.columntable(df)
@time mf = ModelFrame(F, cols, model=LinearModel)
@time mm = ModelMatrix(mf)
@time y = response(mf)
@time linmodel = fit(LinearModel, mm.m, y)
@time regmodel = StatsModels.TableRegressionModel(linmodel, mf, mm);

# %%
typeof(cols)

# %%
typeof(mf)

# %%
typeof(mm) |> x -> (fieldnames(x), fieldtypes(x))

# %%
typeof(y)

# %%
typeof(linmodel)

# %%
typeof(regmodel)

# %%
regmodel.model == linmodel

# %%
regmodel.mf == mf

# %%
regmodel.mm == mm

# %%
@which lm(F, df)

# %%
@which fit(LinearModel, F, df)

# %%
@which fit(LinearModel, mm.m, y)

# %%
@code_warntype term(:y) ~ sum(term(x) for x in x_vars)

# %%
@code_warntype Tables.columntable(df)

# %%
@code_warntype ModelFrame(F, cols, model=LinearModel)

# %%
@code_warntype ModelMatrix(mf)

# %%
@code_warntype response(mf)

# %%
@code_warntype fit(LinearModel, mm.m, y)

# %%
@code_warntype StatsModels.TableRegressionModel(linmodel, mf, mm)

# %%
@code_warntype my_lm(:y, x_vars, df)

# %%

max_stars_repo_path: 0009/fit is very slow for new formulas.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using CPLEX
using Test

function runtests(mod)
    for name in filter(s -> startswith("$(s)", "test_"), names(mod; all = true))
        @testset "$(name)" begin
            getfield(mod, name)()
        end
    end
end

@testset "MathOptInterface Tests" begin
    for file in readdir("MathOptInterface")
        include(joinpath("MathOptInterface", file))
    end
end

@testset "Deprecated functions" begin
    err = ErrorException(CPLEX._DEPRECATED_ERROR_MESSAGE)
    @test_throws err newlongannotation()
    @test_throws err CPLEX.get_status()
    model = CPLEX.Optimizer()
    @test_throws err model.inner
    @test model.lp isa Ptr{Cvoid}
end


text:
"""
    function screen(
        model::MetabolicModel;
        variants::Maybe{Array{V,N}} = nothing,
        analysis,
        args::Maybe{Array{T,N}} = nothing,
        workers = [myid()],
    )::Array where {V<:AbstractVector, T<:Tuple,N}

Take an array of model-modifying function vectors in `variants`, and execute
the function `analysis` on all variants of the `model` specified by `variants`.
The computation is distributed over worker IDs in `workers`. If `args` are
supplied (as an array of the same size as the `variants`), they are forwarded
as arguments to the corresponding analysis function calls.

The array of variants must contain vectors of single-parameter functions, these
are applied to model in order. The functions must *not* modify the model, but
rather return a modified copy. The copy should be made as shallow as possible,
to increase memory efficiency of the process. Variant generators that modify
the argument model in-place will cause unpredictable results. Refer to the
definition of [`screen_variant`](@ref) for details.

The function `analysis` will receive a single argument (the modified model),
together with an expanded tuple of arguments from `args`.

The modification and analysis functions are transferred to `workers` as-is; all
packages required to run them (e.g. the optimization solvers) must be loaded
there. Typically, you want to use the macro `@everywhere using
MyFavoriteSolver` from `Distributed` package for loading the solvers.

# Return value

The results of running `analysis` are collected in to the resulting array, in a
way that preserves the shape of the `variants`, similarly as with `pmap`.

The results of `analysis` function must be serializable, preferably made only
from pure Julia structures, because they may be transferred over the network
between the computation nodes. For that reason, functions that return whole
JuMP models that contain pointers to allocated C structures (such as
[`flux_balance_analysis`](@ref) used with `GLPK` or `Gurobi` otimizers) will
generally not in this context.

# Example
```
function reverse_reaction(i::Int)
    (model::CoreModel) -> begin
        mod = copy(model)
        mod.S[:,i] .*= -1   # this is unrealistic but sufficient for demonstration
        mod
    end
end

m = load_model(CoreModel, "e_coli_core.xml")

screen_variants(m,
           [
               [reverse_reaction(5)],
               [reverse_reaction(3), reverse_reaction(6)]
           ],
           mod -> mod.S[:,3])  # observe the changes in S

screen_variants(m,
    [
        [reverse_reaction(5)],
        [reverse_reaction(3), reverse_reaction(6)]
    ],
    mod -> flux_balance_analysis_vec(mod, GLPK.Optimizer))  # run analysis
"""
function screen(
    model::MetabolicModel;
    variants::Maybe{Array{V,N}} = nothing,
    analysis,
    args::Maybe{Array{T,N}} = nothing,
    workers = [myid()],
)::Array where {V<:AbstractVector,T<:Tuple,N}

    map(fetch, save_at.(workers, :cobrexa_screen_variants_model, Ref(model)))
    map(fetch, save_at.(workers, :cobrexa_screen_variants_analysis_fn, Ref(analysis)))
    map(fetch, get_from.(workers, Ref(:(precache!(cobrexa_screen_variants_model)))))

    if isnothing(variants)
        if isnothing(args)
            throw(
                DomainError(
                    args,
                    "at least one of `variants` and `args` must be non-empty",
                ),
            )
        end
        variants = [[] for _ in args]
    elseif isnothing(args)
        args = [() for _ in variants]
    else
        if size(variants) != size(args)
            throw(
                DomainError(
                    "$(size(variants)) != $(size(args))",
                    "sizes of `variants` and `args` differ",
                ),
            )
        end
    end

    res = dpmap(
        (vars, args)::Tuple -> :($COBREXA.screen_variant(
            cobrexa_screen_variants_model,
            $vars,
            cobrexa_screen_variants_analysis_fn,
            $args,
        )),
        CachingPool(workers),
        zip(variants, args),
    )

    map(fetch, remove_from.(workers, :cobrexa_screen_variants_model))
    map(fetch, remove_from.(workers, :cobrexa_screen_variants_analysis_fn))

    return res
end

"""
    screen_variant(model::MetabolicModel, variant::Vector, analysis, args = ())

Helper function for [`screen`](@ref) that applies all single-argument
functions in `variant` to the `model` (in order from "first" to
"last"), and executes `analysis` on the result.

Can be used to test model variants locally.
"""
function screen_variant(model::MetabolicModel, variant::Vector, analysis, args = ())
    for fn in variant
        model = fn(model)
    end
    analysis(model, args...)
end

"""
    screen_variants(model, variants, analysis; workers=[myid()])

A shortcut for [`screen`](@ref) that only works with model variants.
"""
screen_variants(model, variants, analysis; workers = [myid()]) =
    screen(model; variants = variants, analysis = analysis, workers = workers)


text:
# ======================= SHOWING ====================

Base.show(io::IO, x::PatVar) = print(io, x.name)

function Base.show(io::IO, x::PatLiteral)
    if x.val isa Symbol
        print(io, ":")
    end
    print(io, x.val)
end

Base.show(io::IO, x::PatTypeAssertion) = print(io, x.var, "::", x.type)

Base.show(io::IO, x::PatSplatVar) = print(io, x.name, "...")

Base.show(io::IO, x::PatEquiv) = print(io, x.left, "≡ₙ", x.right)

function Base.show(io::IO, x::PatTerm)
    if x.head == :call
        @assert x.args[1] isa PatLiteral
        print(io, Expr(x.head, x.args[1].val, x.args[2:end]...))
    else
        print(io, Expr(x.head, x.args...))
    end
end

function Base.show(io::IO, x::PatAllTerm)
    n = length(x.args)

    # TODO change me ?
    print(io, "~", x.head)
    print(io, "(")
    for i ∈ 1:n
        @inbounds print(io, x.args[i])
        if i < n
            print(io, ",")
        end
    end
    print(io, ")")
end

# ======================= READING ====================

# Resolve `GlobalRef` instances to literals.
resolve(gr::GlobalRef) = getproperty(gr.mod, gr.name)
resolve(gr) = gr

function Pattern(ex::Expr, mod=@__MODULE__, resolve_fun=false)
    ex = cleanast(ex)
    head = gethead(ex)
    args = getargs(ex)

    n = length(args)
    patargs = Vector{Pattern}(undef, n)
    for i ∈ 1:n
        @inbounds patargs[i] = Pattern(args[i], mod, resolve_fun)
    end

    if head == :call
        if resolve_fun
            f = resolve(GlobalRef(mod, args[1]))
            patargs[1] = PatLiteral(f)
        else
            patargs[1] = PatLiteral(args[1])
        end
    elseif head == :(::)
        v = patargs[1]
        t = patargs[2]
        ty = Union{}
        if t isa PatVar
            ty = getfield(mod, t.name)
        elseif t isa PatLiteral{<:Type}
            ty = t.val
        end
        return PatTypeAssertion(v, ty)
    elseif head == :(...)
        v = patargs[1]
        if v isa PatVar
            return PatSplatVar(v)
        end
    end

    PatTerm(head, patargs)
end

function Pattern(x::Symbol, mod=@__MODULE__, resolve_fun=false)
    PatVar(x)
end

function Pattern(x::QuoteNode, mod=@__MODULE__, resolve_fun=false)
    if x.value isa Symbol
        PatLiteral(x.value)
    else
        PatLiteral(x)
    end
end

# Generic fallback
function Pattern(ex, mod=@__MODULE__, resolve_fun=false)
    if ex isa Expr
        ex = cleanast(ex)
    end
    if isterm(typeof(ex))
        head = gethead(ex)
        args = getargs(ex)

        n = length(args)
        patargs = Vector{Pattern}(undef, n)
        for i ∈ 1:n
            @inbounds patargs[i] = Pattern(args[i], mod, resolve_fun)
        end

        return PatTerm(head, patargs)
    end
    return PatLiteral(ex)
end

function Pattern(p::Pattern, mod=@__MODULE__, resolve_fun=false)
    p
end

macro pat(ex)
    Pattern(ex, __module__, false)
end

macro pat(ex, resolve_fun::Bool)
    Pattern(ex, __module__, resolve_fun)
end


text:
@testset "default options" begin
    method = :CICO_ONE_PASS
    res0 = [get_endpoint(
        [3., 2., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im_dep(x),
        method;
        loss_crit = 9.
    ) for i in 1:3]

    @test isapprox(res0[1].value, 5.0, atol=1e-2)
    @test length(res0[1].profilePoints) > 0
    @test res0[1].status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[1].direction == :right
    @test isapprox(res0[2].value, 2.0+2.0*sqrt(2.), atol=1e-2)
    @test length(res0[2].profilePoints) > 0
    @test res0[2].status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[2].direction == :right
    @test length(res0[3].profilePoints) == 0
    @test res0[3].status == :SCAN_BOUND_REACHED
    @test res0[3].direction == :right
end

@testset ":left" begin
    method = :CICO_ONE_PASS
    res0 = [get_endpoint(
        [3., 2., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im_dep(x),
        method,
        :left;
        loss_crit = 9.,
    ) for i in 1:3]

    @test isapprox(res0[1].value, 1.0, atol=1e-2)
    @test length(res0[1].profilePoints) > 0
    @test res0[1].status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[1].direction == :left
    @test isapprox(res0[2].value, 2.0-2.0*sqrt(2.), atol=1e-2)
    @test length(res0[2].profilePoints) > 0
    @test res0[2].status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[2].direction == :left
    @test length(res0[3].profilePoints) == 0
    @test res0[3].status == :SCAN_BOUND_REACHED
    @test res0[3].direction == :left
end

@testset ":log" begin
    method = :CICO_ONE_PASS
    res0 = [get_endpoint(
        [3., 2., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im_dep(x),
        method,
        :right;
        loss_crit = 9.,
        scale = [:log, :direct, :log]
    ) for i in 1:3]

    @test isapprox(log10(res0[1].value), log10(5.), atol=1e-2)
    @test length(res0[1].profilePoints) > 0
    @test res0[1].status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[1].direction == :right
    @test isapprox(res0[2].value, 2.0+2.0*sqrt(2.), atol=1e-2)
    @test length(res0[2].profilePoints) > 0
    @test res0[2].status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[2].direction == :right
    @test length(res0[3].profilePoints) == 0
    @test res0[3].status == :SCAN_BOUND_REACHED
    @test res0[3].direction == :right
end


text:
## CachedOperator

mutable struct CachedArray{T,N,DM<:AbstractArray{T,N},M<:AbstractArray{T,N}} <: AbstractArray{T,N}
    data::DM
    array::M
    datasize::NTuple{N,Int}
end

const CachedVector{T,DM<:AbstractVector{T},M<:AbstractVector{T}} = CachedArray{T,1,DM,M}
const CachedMatrix{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} = CachedArray{T,2,DM,M}

# CachedArray(data::AbstractArray{T,N}, array::AbstractArray{T,N}, sz::NTuple{N,Int}) where {T,N} =
#     CachedArray{T,N,typeof(data),typeof(array)}(data, array, sz)
CachedArray(data::AbstractArray, array::AbstractArray) = CachedArray(data, array, size(data))

# function CachedArray(::Type{Diagonal}, array::AbstractMatrix{T}) where T
#     axes(array,1) == axes(array,2) || throw(DimensionMismatch("Matrix must be square to cache as diagonal"))
#     CachedArray(Diagonal(Vector{T}(undef, size(array,1))), array)
# end

CachedArray(::Type{Array}, array::AbstractArray{T,N}) where {T,N} =
    CachedArray(Array{T,N}(undef, ntuple(zero,N)), array)

CachedArray(array::AbstractArray{T,N}) where {T,N} =
    CachedArray(similar(array, ntuple(zero,N)), array)

"""
    cache(array::AbstractArray)

Caches the entries of an array.
"""
cache(::Type{MT}, O::AbstractArray) where {MT<:AbstractArray} = CachedArray(MT,O)
cache(A::AbstractArray) = _cache(MemoryLayout(typeof(A)), A)
_cache(_, O::AbstractArray) = CachedArray(O)
_cache(_, O::CachedArray) = CachedArray(copy(O.data), O.array, O.datasize)
_cache(::AbstractStridedLayout, O::AbstractArray) = copy(O)

convert(::Type{AbstractArray{T}}, S::CachedArray{T}) where {T} = S
convert(::Type{AbstractArray{T}}, S::CachedArray) where {T} =
    CachedArray(convert(AbstractArray{T}, S.data), convert(AbstractArray{T}, S.array), S.datasize)

size(A::CachedArray) = size(A.array)
length(A::CachedArray) = length(A.array)

@propagate_inbounds function Base.getindex(B::CachedArray{T,N}, kj::Vararg{Integer,N}) where {T,N}
    @boundscheck checkbounds(B, kj...)
    resizedata!(B, kj...)
    B.data[kj...]
end

@propagate_inbounds function Base.getindex(B::CachedArray{T,1}, k::Integer) where T
    @boundscheck checkbounds(B, k)
    resizedata!(B, k)
    B.data[k]
end

@propagate_inbounds function Base.setindex!(B::CachedArray{T,N}, v, kj::Vararg{Integer,N}) where {T,N}
    @boundscheck checkbounds(B, kj...)
    resizedata!(B,kj...)
    @inbounds B.data[kj...] = v
    v
end

_maximum(ax, I) = maximum(I)
_maximum(ax, ::Colon) = maximum(ax)
function getindex(A::CachedArray, I...)
    @boundscheck checkbounds(A, I...)
    resizedata!(A, _maximum.(axes(A), I)...)
    A.data[I...]
end

getindex(A::CachedVector, ::Colon) = copy(A)
getindex(A::CachedVector, ::Slice) = copy(A)

function getindex(A::CachedVector, I, J...)
    @boundscheck checkbounds(A, I, J...)
    resizedata!(A, _maximum(axes(A,1), I))
    A.data[I]
end

function getindex(A::CachedVector, I::CartesianIndex)
    resizedata!(A, Tuple(I)...)
    A.data[I]
end

function getindex(A::CachedArray, I::CartesianIndex)
    resizedata!(A, Tuple(I)...)
    A.data[I]
end

## Array caching

function resizedata!(B::CachedVector{T,Vector{T}}, n::Integer) where T<:Number
    @boundscheck checkbounds(Bool, B, n) || throw(ArgumentError("Cannot resize beyound size of operator"))

    # increase size of array if necessary
    olddata = B.data
    ν, = B.datasize
    n = max(ν,n)
    if n > length(B.data) # double memory to avoid O(n^2) growing
        B.data = Array{T}(undef, min(2n,length(B.array)))
        B.data[axes(olddata,1)] = olddata
    end

    inds = ν+1:n
    B.data[inds] .= view(B.array,inds)

    B.datasize = (n,)

    B
end

function resizedata!(B::CachedArray{T,N,Array{T,N}},nm::Vararg{Integer,N}) where {T<:Number,N}
    @boundscheck checkbounds(Bool, B, nm...) || throw(ArgumentError("Cannot resize beyound size of operator"))

    # increase size of array if necessary
    olddata = B.data
    νμ = size(olddata)
    nm = max.(νμ,nm)
    if νμ ≠ nm
        B.data = Array{T}(undef, nm)
        B.data[axes(olddata)...] = olddata
    end

    for k in 1:N-1
        inds = tuple(axes(B.data)[1:k-1]...,νμ[k]+1:nm[k],Base.OneTo.(B.datasize[k+1:end])...)
        B.data[inds...] .= view(B.array,inds...)
    end
    let k = N
        inds = tuple(axes(B.data)[1:k-1]...,νμ[k]+1:nm[k])
        B.data[inds...] .= view(B.array,inds...)
    end
    B.datasize = nm

    B
end

function convexunion(a::AbstractVector, b::AbstractVector)
    isempty(a) && return b
    isempty(b) && return a
    min(minimum(a),minimum(b)):max(maximum(a),maximum(b))
end

function colsupport(A::CachedMatrix, i)
    isempty(i) && return 1:0
    minimum(i) ≤ A.datasize[2] ? convexunion(colsupport(A.array, i),colsupport(A.data,i) ∩ Base.OneTo(A.datasize[1])) : colsupport(A.array, i)
end
colsupport(A::CachedVector, i) =
    convexunion(colsupport(A.array, i),colsupport(A.data,i) ∩ Base.OneTo(A.datasize[1]))
function rowsupport(A::CachedMatrix, i)
    isempty(i) && return 1:0
    minimum(i) ≤ A.datasize[1] ? convexunion(rowsupport(A.array, i),rowsupport(A.data,i) ∩ Base.OneTo(A.datasize[2])) : rowsupport(A.array, i)
end

replace_in_print_matrix(A::CachedMatrix, i::Integer, j::Integer, s::AbstractString) =
    i in colsupport(A,j) ? s : replace_with_centered_mark(s)

###
# special for zero cache
###

zero!(A::CachedArray{<:Any,N,<:Any,<:Zeros}) where N = zero!(A.data)
function getindex(A::CachedVector{T,<:AbstractVector,<:AbstractFill{<:Any,1}}, I::AbstractVector) where T
    @boundscheck checkbounds(A, I)
    CachedArray(A.data[I ∩ OneTo(A.datasize[1])], A.array[OneTo(length(I))])
end

###
# MemoryLayout
####

cachedlayout(_, _) = UnknownLayout()
MemoryLayout(C::Type{CachedArray{T,N,DAT,ARR}}) where {T,N,DAT,ARR} = cachedlayout(MemoryLayout(DAT), MemoryLayout(ARR))

#####
# broadcasting
#
# We want broadcasting for numbers with concaenations to pass through
# to take advantage of special implementations of the sub-components
######

BroadcastStyle(::Type{<:CachedArray{<:Any,N}}) where N = LazyArrayStyle{N}()

broadcasted(::LazyArrayStyle, op, A::CachedArray) =
    CachedArray(broadcast(op, paddeddata(A)), broadcast(op, A.array))

broadcasted(::LazyArrayStyle, op, A::CachedArray, c::Number) =
    CachedArray(broadcast(op, paddeddata(A), c), broadcast(op, A.array, c))
broadcasted(::LazyArrayStyle, op, c::Number, A::CachedArray) =
CachedArray(broadcast(op, c, paddeddata(A)), broadcast(op, c, A.array))
broadcasted(::LazyArrayStyle, op, A::CachedArray, c::Ref) =
    CachedArray(broadcast(op, paddeddata(A), c), broadcast(op, A.array, c))
broadcasted(::LazyArrayStyle, op, c::Ref, A::CachedArray) =
    CachedArray(broadcast(op, c, paddeddata(A)), broadcast(op, c, A.array))

function cache_broadcast(op, A::CachedVector, B)
    dat = paddeddata(A)
    n = length(dat)
    m = length(B)
    CachedArray(broadcast(op, dat, view(B,1:n)), broadcast(op, A.array, B))
end

function cache_broadcast(op, A, B::CachedVector)
    dat = paddeddata(B)
    n = length(dat)
    m = length(A)
    CachedArray(broadcast(op, view(A,1:n), dat), broadcast(op, A, B.array))
end

broadcasted(::LazyArrayStyle, op, A::CachedVector, B::AbstractVector) = cache_broadcast(op, A, B)
broadcasted(::LazyArrayStyle, op, A::AbstractVector, B::CachedVector) = cache_broadcast(op, A, B)

function broadcasted(::LazyArrayStyle, op, A::CachedVector, B::CachedVector)
    n = max(A.datasize[1],B.datasize[1])
    resizedata!(A,n)
    resizedata!(B,n)
    Adat = view(paddeddata(A),1:n)
    Bdat = view(paddeddata(B),1:n)
    CachedArray(broadcast(op, Adat, Bdat), broadcast(op, A.array, B.array))
end

###
# norm
###

norm1(a::CachedVector) = norm(paddeddata(a),1) + norm(@view(a.array[a.datasize[1]+1:end]),1)
norm2(a::CachedVector) = sqrt(norm(paddeddata(a),2)^2 + norm(@view(a.array[a.datasize[1]+1:end]),2)^2)
normInf(a::CachedVector) = max(norm(paddeddata(a),Inf), norm(@view(a.array[a.datasize[1]+1:end]),Inf))
normp(a::CachedVector, p) = (norm(paddeddata(a),2)^p + norm(@view(a.array[a.datasize[1]+1:end]),2)^p)^inv(p)


text:
function ptcKrylovTestB(n = 63)
    delta0 = .01
    pout1 = ptciBeam(;lsolver="bicgstab");
    pout2 = ptciBeam(n, delta0, false; lsolver="bicgstab");
    pout3 = ptciBeam(n, delta0, false, "left"; lsolver="bicgstab");
    sol1 = pout1.solution
    sol2 = pout2.solution
    sol3 = pout3.solution
#
# sol3 is the wrong stable branch. Left preconditioning bites you!
#
solpass1a = (norm(sol1 - sol2, Inf) < 1.e-9)
solpass1b = (norm(sol1 - sol3, Inf) < 1.e-9)
solpass1 = solpass1a && solpass1b
solpass1 || println("solpass1 fails for ptcsoli-bicgstab")
histpass = ( length(pout1.history) == 25 )
histpass || println("histpass fails for ptcsoli-bicgstab")
solpass2 = ( abs(norm(sol1,Inf) - 2.191) < 1.e-3)
solpass2 || println("solpass2 fails for ptcsoli-bicgstab")
ptcipass = solpass1 && histpass && solpass2

end


text:
import Base: +, -, *, /, ==, isless, rem, convert, zero, one, iseven
import Base: show

"""
    Phase
The type supports manipulating phases as expressions.
`Phase(x)` represents the number `x⋅π`.
"""
struct Phase
    ex
    type
end

Phase(p::T) where {T} = Phase(p, T)

Phase(p::Phase) = p

function show(io::IO, p::Phase)
    if p.ex isa Number
        print(io, "$(p.ex)⋅π")
    else
        print(io, "Phase(($(p.ex))::$(p.type))")
    end
end

function +(p1::Phase, p2::Phase)
    T1 = p1.type
    T2 = p2.type
    if p1.ex isa Number && p2.ex isa Number
        return Phase(p1.ex + p2.ex)
    end

    T = Base.promote_op(+, T1, T2)
    return Phase(Expr(:call, :+, p1.ex, p2.ex), T)
end
+(p1::Phase, p2::Number) = p1 + Phase(p2)
+(p1::Number, p2::Phase) = Phase(p1) + p2

function -(p1::Phase, p2::Phase)
    T1 = p1.type
    T2 = p2.type
    if p1.ex isa Number && p2.ex isa Number
        return Phase(p1.ex - p2.ex)
    end

    T = Base.promote_op(-, T1, T2)
    return Phase(Expr(:call, :-, p1.ex, p2.ex), T)
end
-(p1::Phase, p2::Number) = p1 - Phase(p2)
-(p1::Number, p2::Phase) = Phase(p1) - p2

function *(p1::Phase, p2::Phase)
    T1 = p1.type
    T2 = p2.type
    if p1.ex isa Number && p2.ex isa Number
        return Phase(p1.ex * p2.ex)
    end

    T = Base.promote_op(*, T1, T2)
    return Phase(Expr(:call, :*, p1.ex, p2.ex), T)
end
*(p1::Phase, p2::Number) = p1 * Phase(p2)
*(p1::Number, p2::Phase) = Phase(p1) * p2

function /(p1::Phase, p2::Phase)
    T1 = p1.type
    T2 = p2.type
    if p1.ex isa Number && p2.ex isa Number
        return Phase(p1.ex / p2.ex)
    end

    T = Base.promote_op(/, T1, T2)
    return Phase(Expr(:call, :/, p1.ex, p2.ex), T)
end
/(p1::Phase, p2::Number) = p1 / Phase(p2)
/(p1::Number, p2::Phase) = Phase(p1) / p2

function -(p::Phase)
    T0 = p.type
    if p.ex isa Number
        return Phase(-p.ex)
    end

    T = Base.promote_op(-, T0)
    return Phase(Expr(:call, :-, p.ex), T)
end

==(p1::Phase, p2::Phase) = p1.ex == p2.ex
==(p1::Phase, p2::Number) = (p1 == Phase(p2))
==(p1::Number, p2::Phase) = (Phase(p1) == p2)

isless(p1::Phase, p2::Number) = (p1.ex isa Number) && p1.ex < p2
function rem(p::Phase, d::Number)
    if p.ex isa Number
        return Phase(rem(p.ex, d))
    end
    return p
end

convert(::Type{Phase}, p) = Phase(p)
convert(::Type{Phase}, p::Phase) = p

zero(::Phase) = Phase(0//1)
zero(::Type{Phase}) = Phase(0//1)
one(::Phase) = Phase(1//1)
one(::Type{Phase}) = Phase(1//1)

iseven(p::Phase) = (p.ex isa Number) && (-1)^p.ex > 0


text:
export checkDerivative

"""
checkDerivative(f,x0;out,tol,nSuccess)

Check if function f provides correct derivatives around point x0 using a Taylor expansion.

Inputs:

	f::Function  - fx,dfx= f(x) provides function value and derivative
	x::Array     - vector around which to test the derivatives
	out::Bool    - flag for activating printing (default: true)
	tol::Number  - tolerance to judge if error in linearization decreases quadratically (default: 1.9)
	nSuccess:Int - number of quadratically-decreasing steps needed to decide that derivative is correct

Outputs:
	pass, Error, Order

"""
function checkDerivative(f::Function,x0;out::Bool=true,tol::Number=1.9,nSuccess::Int=3)
	if out
		println(@sprintf("%9s\t%9s\t%9s\t%9s\t%9s\t%5s","h","E0","E1","O1","O2","OK?"))
	end
	v = randn(size(x0))
	if eltype(x0) <: Complex
		v += 1im*randn(size(v))
	end
	f0,dvf = f(x0,v)
	dvf   = real(dvf)
	Error = zeros(10,2)
	Order = zeros(10,2)
	Success = zeros(10)
	for j=1:10
		ft = f(x0+10.0^(-j)*v)            # function value
		Error[j,1] = norm(f0-ft)/norm(f0)           # Error TaylorPoly 0
		Error[j,2] = norm(f0 .+10.0^(-j)*dvf .- ft)/norm(f0) # Error TaylorPoly 1
		if j>1
			Order[j,:] = log10.(Error[j-1,:]./Error[j,:]);
		end
		if (Order[j,2]>tol) || (Error[j,1]/Error[j,2] > 100); Success[j]=1; end
		if out
			println(@sprintf("%1.3e\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t%5d",
							10.0^(-j), Error[j,1],Error[j,2], Order[j,1],Order[j,2],Success[j]))
		end
	end
	pass = sum(Success) > nSuccess
	return  pass,Error,Order
end


text:
using MyResults;

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

include("res_repo_tests/res_repo_tests.jl");


text:
function default_min_tol(est; n_min=20)
    n_min = min(n_min, length(est.eps))
    est.eps[n_min]
end

"""
    plot_abc(est; min_tol=default, gather=true, labels=nothing, xlabel="\\epsilon")

Plot ABC post-correction estimators with related 95% confidence intervals.

# Arguments
- `est`: Output of `abc_postprocess`.
- `min_tol`: Minimum tolerance. By default, 20th smallest sample tolerance.
- `gather`: Whether the individual plots are gathered to one.
- `labels`: The labels of the estimates.
- `xlabel`: The label of the x-axis (tolerance).
"""
function plot_abc(est; min_tol=default_min_tol(est), gather=true, labels=nothing,
    xlabel="\\epsilon")
    lims(x) = (minimum(x), maximum(x))
    d, n = size(est.E)
    v = est.eps .>= min_tol
    eps_ = est.eps[v]
    if length(eps_) == 0
        warning("All distances less than min_tol")
        return nothing
    end
    p = Vector{Plots.Plot}(undef, 0)
    for i in 1:d
        E_ = est.E[i,v]; L_ = est.ci_L[i,v];  U_ = est.ci_U[i,v]
        push!(p, Plots.plot(eps_, E_, ribbon=[E_-L_, U_-E_], color=:black,
                          fillalpha=0.3, legend=false, xlim=lims(eps_),
                          framestyle=:axes))
        if !(labels isa Nothing)
            Plots.plot!(p[end], ylabel=labels[i])
        end
    end
    if gather
        Plots.plot!(p[end], xlabel=xlabel)
        Plots.plot(p..., layout=(d,1))
    else
        p
    end
end


text:
# This file contains functions to calculate an analytical solution of Riemann problem.
#
# Reference: Toro (2013)

function P_star()
  # initial guess
  P0 = 0.5 * (P_left + P_right)

  # Newton-Raphson iterative procedure
  CHA = 1
  while CHA > TOL
    P1 = P0 - f(P0) / dfdP(P0)
    CHA = abs(P1 - P0) / 0.5 / (P1 + P0)
    P0 = P1
    @show P0
  end
  return P0
end

function v_star(P)
  return 0.5 * (v_left + v_right) + 0.5(f_right(P) - f_left(P))
end

function region(x, t, v, P, x_center = 0)
  # 1: left
  # 2: right
  # 3: left star (shock)
  # 4: right star (shock)
  # 5: left star (rarefaction)
  # 6: right star (rarefaction)
  # 7: left rarefaction fan
  # 8: right rarefaction fan

  S = (x - x_center) / t
  if S < v # left from contact discontinuity
    if P > P_left # shock
      if S < S_shock_left(P) # outside shock
        return 1 #"left"
      else # inside shock
        return 3 #"left star"
      end
    else # rarefaction
      if S < S_head_left() # outside
        return 1 #"left"
      elseif S < S_tail_left(P) # inside rarefaction fan
        return 7 #"left rarefaction fan"
      else
        return 5 #"left star"
      end
    end
  else # right from contact discontinuity
    if P > P_right
      if S > S_shock_right(P)
        return 2 #"right"
      else
        return 4 #"right star"
      end
    else
      if S > S_head_right()
        return 2 #"right"
      elseif S > S_tail_right(P)
        return 8 #"right rarefaction fan"
      else
        return 6 #"right star"
      end
    end
  end
end

function rho_riemann(x, t, v, P, x_center = 0)
  reg = region(x, t, v, P, x_center)
  if reg == 1 #"left"
    return rho_left
  elseif reg == 2 #"right"
    return rho_right
  elseif reg == 3 # left star (shock)
    return rho_left * (P / P_left + (gamma - 1) / (gamma + 1)) /
           ((gamma - 1) / (gamma + 1) * P / P_left + 1)
  elseif reg == 4 # right star (shock)
    return rho_right * (P / P_right + (gamma - 1) / (gamma + 1)) /
           ((gamma - 1) / (gamma + 1) * P / P_right + 1)
  elseif reg == 5 #"left star (rarefaction)"
    return rho_left * (P / P_left)^(1 / gamma)
  elseif reg == 6 #"right star (rarefaction"
    return rho_right * (P / P_right)^(1 / gamma)
  elseif reg == 7 #"left rarefaction fan"
    return rho_left *
           (
      2 / (gamma + 1) + (gamma - 1) / (gamma + 1) / aL() * (v_left - x / t)
    )^(2 / (gamma - 1))
  elseif reg == 8 #"right rarefaction fan"
    return rho_right *
           (
      2 / (gamma + 1) - (gamma - 1) / (gamma + 1) / aR() * (v_right - x / t)
    )^(2 / (gamma - 1))
  else
    @assert 0
  end
end

function v_riemann(x, t, v, P, x_center = 0)
  reg = region(x, t, v, P, x_center)
  if reg == 1 #"left"
    return v_left
  elseif reg == 2 #"right"
    return v_right
  elseif reg == 3 || reg == 4 || reg == 5 || reg == 6 #"left star" || "right star"
    return v
  elseif reg == 7 #"left rarefaction fan"
    return 2 / (gamma + 1) * (aL() + (gamma - 1) / 2 * v_left + x / t)
  elseif reg == 8 #"right rarefaction fan"
    return 2 / (gamma + 1) * (-aR() + (gamma - 1) / 2 * v_right + x / t)
  else
    @assert 0
  end
end

function P_riemann(x, t, v, P, x_center = 0)
  reg = region(x, t, v, P, x_center)
  if reg == 1 #"left"
    return P_left
  elseif reg == 2 #"right"
    return P_right
  elseif reg == 3 || reg == 4 || reg == 5 || reg == 6 #"left star" || "right star"
    return P
  elseif reg == 7 #"left rarefaction fan"
    return P_left *
           (
      2 / (gamma + 1) + (gamma - 1) / (gamma + 1) / aL() * (v_left - x / t)
    )^(2 * gamma / (gamma - 1))
  elseif reg == 8 #"right rarefaction fan"
    return P_right *
           (
      2 / (gamma + 1) - (gamma - 1) / (gamma + 1) / aR() * (v_right - x / t)
    )^(2 * gamma / (gamma - 1))
  else
    @assert 0
  end
end

function f_left(P)
  if P > P_left # shock
    f = (P - P_left) * sqrt(AL() / (P + BL()))
  else # rarefaction
    f = 2 * aL() / (gamma - 1) * ((P / P_left)^((gamma - 1) / 2 / gamma) - 1)
  end
  return f
end

function f_right(P)
  if P > P_right # shock
    f = (P - P_right) * sqrt(AR() / (P + BR()))
  else # rarefaction
    f = 2 * aR() / (gamma - 1) * ((P / P_right)^((gamma - 1) / 2 / gamma) - 1)
  end
  return f
end

function dfdP_left(P)
  if P > P_left
    dfdP = sqrt(AL() / (BL() + P)) * (1 - (P - P_left) / 2 / (BL() + P))
  else
    dfdP = 1 / rho_left / aL() * (P / P_left)^(-(gamma + 1) / 2 / gamma)
  end
  return dfdP
end

function dfdP_right(P)
  if P > P_right
    dfdP = sqrt(AR() / (BR() + P)) * (1 - (P - P_right) / 2 / (BR() + P))
  else
    dfdP = 1 / rho_right / aR() * (P / P_right)^(-(gamma + 1) / 2 / gamma)
  end
  return dfdP
end

function AL()
  return 2 / (gamma + 1) / rho_left
end
function AR()
  return 2 / (gamma + 1) / rho_right
end
function BL()
  return (gamma - 1) / (gamma + 1) * P_left
end
function BR()
  return (gamma - 1) / (gamma + 1) * P_right
end
function aL()
  return sqrt(gamma * P_left / rho_left)
end
function aR()
  return sqrt(gamma * P_right / rho_right)
end

function f(P)
  return f_left(P) + f_right(P) + v_right - v_left
end

function dfdP(P)
  return dfdP_left(P) + dfdP_right(P)
end

function S_shock_left(P)
  return v_left -
         aL() * ((gamma + 1) / 2 / gamma * P / P_left + (gamma - 1) / 2 / gamma)^(1 / 2)
end

function S_shock_right(P)
  return v_right +
         aR() * ((gamma + 1) / 2 / gamma * P / P_right + (gamma - 1) / 2 / gamma)^(1 / 2)
end

function S_head_left()
  return v_left - aL()
end
function S_tail_left(P)
  rho_s = rho_left * (P / P_left)^(1 / gamma)
  v_s = v_star(P)
  aL_s = sqrt(gamma * P / rho_s)
  return v_s - aL_s
end
function S_head_right()
  return v_right + aR()
end
function S_tail_right(P)
  rho_s = rho_right * (P / P_right)^(1 / gamma)
  v_s = v_star(P)
  aR_s = sqrt(gamma * P / rho_s)
  return v_s + aR_s
end


text:
# Autogenerated wrapper script for hsakmt_roct_jll for x86_64-linux-gnu
export libhsakmt

using NUMA_jll
JLLWrappers.@generate_wrapper_header("hsakmt_roct")
JLLWrappers.@declare_library_product(libhsakmt, "libhsakmt.so.1")
function __init__()
    JLLWrappers.@generate_init_header(NUMA_jll)
    JLLWrappers.@init_library_product(
        libhsakmt,
        "lib/libhsakmt.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
## Test fit with only one hidden state, fit! should equal MLE
hmm = HMM(1,Normal())

# generate some random synthetic data
μ = randn()*10
σ = 1 + rand()*5
o = μ + randn(100)*σ

fit!(hmm, o; tol=NaN)

# use mle, not un-biased estimators
μ_true = mean(o)
σ_true = sqrt(sum((o - mean(o)).^2) / length(o))

@test(round(hmm.B[1].μ,10) == round(μ_true,10))
@test(round(hmm.B[1].σ,10) == round(σ_true,10))

## Test fit on easy synthetic data with two hidden states
μ1,μ2 = -20,15
σ1,σ2 = 3,5

A = [ 0.9 0.1 ;
      0.8 0.2 ]

B = (Distribution)[ Normal(μ1,σ1) , Normal(μ2,σ2) ]

hmm_true = HMM(A,B)
s,o = generate(hmm_true,10_000)

hmm = HMM(2,Normal())
fit!(hmm,o)

@test( (abs(hmm.B[1].μ - μ1) < 0.5 && abs(hmm.B[1].σ - σ1) < 0.5) ||
       (abs(hmm.B[2].μ - μ1) < 0.5 && abs(hmm.B[2].σ - σ1) < 0.5) )
@test( (abs(hmm.B[1].μ - μ2) < 0.5 && abs(hmm.B[1].σ - σ2) < 0.5) ||
       (abs(hmm.B[2].μ - μ2) < 0.5 && abs(hmm.B[2].σ - σ2) < 0.5) )


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Simulate_PendulumDAE

Simulate PendulumDAE model (= index-3 model transformed to index-1 model)
"""
module Simulate_PendulumDAE

include(joinpath("models", "PendulumDAE.jl"))
import .PendulumDAE

import ModiaMath

model  = PendulumDAE.Model()
result = ModiaMath.simulate!(model, stopTime=2.0, log=true)

ModiaMath.plot(result, [(:x, :y), (:vx, :vy), :lambda], heading="Simulate_PendulumDAE.jl (index3 reduced to index1)")

end


text:
using MultipathChannel
using Test
using Statistics

@testset "MultipathChannel.jl" begin

    @testset "Average power after fading" begin
        maxdoppler = 0
        samplingrate = 15.36e6
        channelmodel = "TDLC"
        rmsdelay = 300e-9

        t = 0:0.05:50
        sinewave = sin.(2π*t)

        numiters = 1000
        power = Vector{Float64}(undef, numiters)
        for i = 1:numiters
            chan = initChan(maxdoppler, samplingrate, channelmodel, rmsdelay=rmsdelay)
            response = applyChan(sinewave, chan)
            power[i] = sum((abs.(response)).^2)/length(sinewave)
        end
        @test abs(mean(power) - mean(sum((abs.(sinewave)).^2)/length(sinewave))) < 0.05
    end

    @testset "Expected value of Rayleigh distribution" begin
        maxdoppler = 0
        samplingrate = 15.36e6
        channelmodel = "rayleigh"

        numiters = 1000
        amp = Vector{Complex{Float64}}(undef, numiters)
        for i = 1:numiters
            chan = initChan(maxdoppler, samplingrate, channelmodel)
            amp[i] = applyChan([1], chan, initialtime=rand()*10000)[1]
        end
        σ = 1/√2
        expected = σ*√(π/2)
        @test abs(mean(abs.(amp)) - expected) < 0.05
    end

end


text:
module openbuildnet

export
# Basics
OBNUpdateMask, OBNSimTimeType,

# Node interface
OBNNode, delete,
on_block_output, on_block_state, on_init, on_term, on_restart,
run, runsafe, process_port_events, stop, schedule,
isstopped, isrunning, iserror,
sim_time, sim_time_ticks, timeunit, wallclock_time, to_ticks,

# Port interface
OBNPort, OBNInputAbstract, OBNOutputAbstract, OBNInput, OBNOutput,
create_input, create_output,
get, set, sendsync, pending, connectfrom, portinfo,
on_receive,

# Misc
lastErrorMessage, max_blockid,
update_mask

include("obndefs.jl")
include("obnextapi.jl")
include("obnnode.jl")
include("obnport.jl")

end # module


text:
using JLD2, LinearAlgebra, Oceananigans, Printf, Statistics
record_interaction = false
include(pwd() * "/scripts/vizinanigans.jl")
include(pwd() * "/scripts/states.jl")
include(pwd() * "/scripts/compare.jl")

include(pwd() * "/analysis_scripts/" * "post_analysis.jl") # Gradients etc. here

files = [pwd() * "/Channel_1_checkpoint_iteration404905.jld2",
         pwd() * "/Channel_2_checkpoint_iteration728850.jld2",
         pwd() * "/Channel_3_checkpoint_iteration1087211.jld2",
         pwd() * "/Channel_4_checkpoint_iteration6160221.jld2",
         pwd() * "/Channel_8_checkpoint_iteration3164301.jld2",
         pwd() * "/Channel_16_checkpoint_iteration6317902.jld2",
         pwd() * "/Channel_24_checkpoint_iteration8612852.jld2",
         pwd() * "/Channel_32_checkpoint_iteration1272125.jld2"
]

weak_files = [
        pwd() * "/Weak_Channel_1_checkpoint_iteration370489.jld2",
        pwd() * "/Weak_Channel_3_checkpoint_iteration1087211.jld2"
]

new_files = [
    pwd() * "/Channel_1_checkpoint_iteration370489.jld2",
    pwd() * "/Channel_3_checkpoint_iteration1075121.jld2",
    pwd() * "/Channel_16_checkpoint_iteration1003423.jld2",
]

ridge_files = [
    pwd() * "/Ridge_8_checkpoint_iteration939858.jld2",
    pwd() * "/Ridge_8_checkpoint_iteration6078160.jld2",
    pwd() * "/Ridge_16_checkpoint_iteration1051200.jld2",
]

relaxation_files = [
    pwd() * "/Relaxation_Channel_16_checkpoint_iteration1445551.jld2"
]

abernathy_files = [
    pwd() * "/Abernathy_16_checkpoint_iteration1158859.jld2"
    pwd() * "/Abernathy_16_checkpoint_iteration286890.jld2"
    pwd() * "/Abernathy_16_τ0.20_Q1_checkpoint_iteration2866921.jld2"
]
##
# http://juliaplots.org/MakieReferenceImages/gallery/index.html
filename = abernathy_files[end]
states, statenames, units = grabstates(filename)
scene = volumeslice(states, statenames = statenames, aspect = (1, 2, 32/192),
                    statistics = true, units = units, statlabelsize = (15, 15) );

# display(scene)
## save interaction
seconds = 20
fps = 10
frames = round(Int, fps * seconds )
if record_interaction
GLMakie.record(scene, pwd() * "/ridge.mp4"; framerate = fps) do io
    for i = 1:frames
        sleep(1/fps)
        recordframe!(io)
    end
end
end
##
filename = new_files[1]
states, statenames, units1 = grabstates(filename)
title = "" * grabtitle(filename)

filename2 = new_files[end]
states2, statenames2, units2 = grabstates(filename2)
title2 = " "  * grabtitle(filename2)

scene = visualize(states, states2, statenames = statenames, statenames2 = statenames2, aspect = (1,1, 32/192), statistics = true, title = title, title2 = title2, units1 = units1, units2 = units2)
##
# instantaneous zonal
# filename = pwd() * "/Weno_20_checkpoint_iteration21030963.jld2" # new_files[1]
states, statenames, units1 = grabstates(filename)
function zonalmean(state)
    return mean(state, dims = 1)[1,:,:]
end
meanstates = zonalmean.(states)
statenames[3]
v̅ =  zonalmean(states[3])
b̅ =  zonalmean(states[5])
vb = zonalmean(states[3] .* states[5])
visualize([v̅, b̅, vb, vb - v̅ .* b̅ ], statenames = ["v̄", "b̄", "avg(vb)",  " v'b' "], title = "Zonal Average")
scene = visualize(meanstates, statenames = statenames, title = "Zonal Average ")
##
filename = new_files[end]
states, statenames, units1 = grabstates(filename)
meanstates2 = zonalmean.(states)
coarsegrainedb = coarsey(meanstates2[5], 12)
coarseb = meanstates[5]
scene = visualize([coarseb, coarsegrainedb, (coarsegrainedb-coarseb) ./ norm(coarseb)], statenames = ["b", "grained b", "difference"], title = "Coarse ")
##
seconds = 20
fps = 30
frames = round(Int, fps * seconds )
if record_interaction
record(scene, pwd() * "/test.mp4"; framerate = fps) do io
    for i = 1:frames
        sleep(1/fps)
        recordframe!(io)
    end
end
end

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
using Test
prec = 12
@testset "soilstorage" begin
    for i in 1:10000
        #parameters that don't change
        beta = rand(0.5:0.1:4)[1]
        Ce =  rand(0.4:0.01:0.8)[1]
        #Percolationcapacity = rand(0.1:0.1:4)[1]
        Ratio_Pref = rand(0.1:0.1:0.9)[1]
        Soilstoragecapacity = rand(50.:5.:400.)[1]

        #parameters that change
        Effective_Precipitation = rand(0.1:0.5:15)[1]
        Potential_Evaporation = rand(0.1:0.2:5)[1]
        Interception_Evaporation = rand(0:0.1:Potential_Evaporation * 0.5)[1]
        Storage = rand(1:5:Soilstoragecapacity)[1]

        # if no precipitation and storage zero, all outputs are zero
        Storage = 0.
        Effective_Precipitation = 0.
        Overlandflow, Preferentialflow, Soil_Evaporation, Soilstorage = soilstorage(Effective_Precipitation, Interception_Evaporation, Potential_Evaporation, Storage, beta, Ce, Ratio_Pref, Soilstoragecapacity)
        @test Overlandflow == 0
        #@test Percolationflow == 0
        @test Preferentialflow == 0
        @test Soil_Evaporation == 0
        @test Soilstorage == 0

        # if precipitation is zero
        Storage = rand(1:5:Soilstoragecapacity)[1]
        Effective_Precipitation = 0.
        Overlandflow, Preferentialflow, Soil_Evaporation, Soilstorage = soilstorage(Effective_Precipitation, Interception_Evaporation, Potential_Evaporation, Storage, beta, Ce, Ratio_Pref, Soilstoragecapacity)
        @test Overlandflow == 0
        @test Preferentialflow == 0
        @test Soil_Evaporation == min(Storage,(Potential_Evaporation - Interception_Evaporation) * min((Storage)/(Soilstoragecapacity * Ce),1))
        #@test Percolationflow == (Storage - Soil_Evaporation)/Soilstoragecapacity * Percolationcapacity
        @test Soilstorage == Storage - Soil_Evaporation

        # if precipitation is not zero
        Storage = rand(1:5:Soilstoragecapacity)[1]
        Effective_Precipitation = rand(0.1:0.5:15)[1]
        Overlandflow, Preferentialflow, Soil_Evaporation, Soilstorage = soilstorage(Effective_Precipitation, Interception_Evaporation, Potential_Evaporation, Storage, beta, Ce, Ratio_Pref, Soilstoragecapacity)
        DischargeSoil = min((1 - (1 - (1 - Storage/Soilstoragecapacity)^beta)) * Effective_Precipitation, Soilstoragecapacity - Storage)
        @test Overlandflow == (Effective_Precipitation - DischargeSoil) * Ratio_Pref
        # next tests sometimes has roundin errors!!!
        @test round(Preferentialflow, digits=prec) == round(Effective_Precipitation - DischargeSoil - Overlandflow, digits=prec)
        Storage_New = Storage + DischargeSoil
        @test Soil_Evaporation == min(Storage_New,(Potential_Evaporation - Interception_Evaporation) * min((Storage_New)/(Soilstoragecapacity * Ce),1))
        #@test Percolationflow == (Storage_New - Soil_Evaporation)/Soilstoragecapacity * Percolationcapacity
        @test Soilstorage == Storage_New - Soil_Evaporation

        # if Storage = Capacity no water flows into soil
        Storage = Soilstoragecapacity
        Effective_Precipitation = rand(0.1:0.5:15)[1]
        Overlandflow, Preferentialflow, Soil_Evaporation, Soilstorage = soilstorage(Effective_Precipitation, Interception_Evaporation, Potential_Evaporation, Storage, beta, Ce, Ratio_Pref, Soilstoragecapacity)
        #DischargeSoil = (1 - (1 - (1 - Storage/Soilstoragecapacity)^beta)) * Effective_Precipitation
        @test round(Overlandflow + Preferentialflow, digits=prec) == round(Effective_Precipitation, digits=prec)
        @test Soil_Evaporation == min(Storage,(Potential_Evaporation - Interception_Evaporation) * min((Storage)/(Soilstoragecapacity * Ce),1))
        #@test Percolationflow == (Storage - Soil_Evaporation)/Soilstoragecapacity * Percolationcapacity
        @test Soilstorage == Storage - Soil_Evaporation

    end
end

using Test

# @testset "soilstorage_rip" begin
#     for i in 1:1000
#         #parameters that don't change
#         beta = rand(0.5:0.1:4)[1]
#         Ce =  rand(0.4:0.01:0.8)[1]
#         Drainagecapacity = rand(0.1:0.1:6)[1]
#         #Ratio_Pref = rand(0.1:0.1:0.9)[1]
#         Soilstoragecapacity = rand(50.:5.:400.)[1]
#
#         #parameters that change
#         Effective_Precipitation = rand(0.1:0.5:15)[1]
#         Potential_Evaporation = rand(0.1:0.2:10)[1]
#         Interception_Evaporation = rand(0:0.1:Potential_Evaporation * 0.5)[1]
#         Riparian_Discharge = rand(0.1:0.2:10)[1]
#         Storage = rand(1:5:Soilstoragecapacity)[1]
#
#         # if no precipitation and storage zero, all outputs are zero
#         Storage = 0.
#         Effective_Precipitation = 0.
#         Riparian_Discharge = 0.
#         Overlandflow, Soil_Evaporation, Soilstorage = ripariansoilstorage(Effective_Precipitation, Interception_Evaporation, Potential_Evaporation, Riparian_Discharge, Storage, beta, Ce, Drainagecapacity, Soilstoragecapacity)
#         @test Overlandflow == 0
#         @test Soil_Evaporation == 0
#         @test Soilstorage == 0
#
#         # if precipitation is zero
#         Storage = rand(1:5:Soilstoragecapacity)[1]
#         Riparian_Discharge = rand(0.1:0.2:10)[1]
#         Effective_Precipitation = 0.
#         Overlandflow, Soil_Evaporation, Soilstorage = ripariansoilstorage(Effective_Precipitation, Interception_Evaporation, Potential_Evaporation, Riparian_Discharge, Storage, beta, Ce, Drainagecapacity, Soilstoragecapacity)
#         # overlandflow occurs due to Riparian Discharge entering the soil
#         Ratio_Soil = (1 - (1 - Storage/Soilstoragecapacity)^beta)
#         DischargeSoil = min((1 - Ratio_Soil) * (Riparian_Discharge), Soilstoragecapacity - Storage)
#         # overland flow also includes the fast drainage!!
#         Storage_New = Storage + DischargeSoil
#         @test Soil_Evaporation == min(Storage_New,(Potential_Evaporation - Interception_Evaporation) * min((Storage_New)/(Soilstoragecapacity * Ce),1))
#         Fastdrainage = ((Storage_New - Soil_Evaporation) / Soilstoragecapacity * Drainagecapacity)
#         @test Soilstorage == Storage + DischargeSoil - Soil_Evaporation - Fastdrainage
#         @test Overlandflow == Riparian_Discharge - DischargeSoil + Fastdrainage
#         # if precipitation is not zero
#         Storage = rand(1:5:Soilstoragecapacity)[1]
#         Effective_Precipitation = rand(0.1:0.5:15)[1]
#         Overlandflow, Soil_Evaporation, Soilstorage = ripariansoilstorage(Effective_Precipitation, Interception_Evaporation, Potential_Evaporation, Riparian_Discharge, Storage, beta, Ce, Drainagecapacity, Soilstoragecapacity)
#         Ratio_Soil = (1 - (1 - Storage/Soilstoragecapacity)^beta)
#         DischargeSoil = min((1 - Ratio_Soil) * (Effective_Precipitation + Riparian_Discharge), Soilstoragecapacity - Storage)
#         Storage_New = Storage + DischargeSoil
#
#         @test Soil_Evaporation == min(Storage_New,(Potential_Evaporation - Interception_Evaporation) * min((Storage_New)/(Soilstoragecapacity * Ce),1))
#         Fastdrainage = ((Storage_New - Soil_Evaporation) / Soilstoragecapacity * Drainagecapacity)
#         @test Soilstorage == Storage_New - Soil_Evaporation - Fastdrainage
#         @test Overlandflow == Riparian_Discharge + Effective_Precipitation - DischargeSoil + Fastdrainage
#
#         # if Storage = Capacity no water flows into soil
#         Storage = Soilstoragecapacity
#         Effective_Precipitation = rand(0.1:0.5:15)[1]
#         Overlandflow, Soil_Evaporation, Soilstorage = ripariansoilstorage(Effective_Precipitation, Interception_Evaporation, Potential_Evaporation, Riparian_Discharge, Storage, beta, Ce, Drainagecapacity, Soilstoragecapacity)
#         @test Soil_Evaporation == min(Storage,(Potential_Evaporation - Interception_Evaporation) * min((Storage)/(Soilstoragecapacity * Ce),1))
#         Fastdrainage = (Storage - Soil_Evaporation)/Soilstoragecapacity * Drainagecapacity
#         @test Soilstorage == Storage - Soil_Evaporation - Fastdrainage
#         @test Overlandflow == Effective_Precipitation + Riparian_Discharge + Fastdrainage
#
#     end
# end
#
# @testset "soilstorage_run" begin
#     for i in 1:100
#         #parameters that don't change
#         beta = rand(0.5:0.1:4)[1]
#         Ce =  rand(0.4:0.01:0.8)[1]
#         #Percolationcapacity = rand(0.1:0.1:4)[1]
#         Ratio_Pref = rand(0.1:0.1:0.9)[1]
#         Soilstoragecapacity = rand(50.:5:400.)[1]
#
#         #parameters that change
#         #Effective_Precipitation = rand(0.1:0.5:15)[1]
#         #Potential_Evaporation = rand(0.1:0.2:10)[1]
#         #Storage = rand(1:5:Soilstoragecapacity)[1]
#
#         # if no precipitation and storage zero, all outputs are zero
#         # test if soilstorage always lower than capacity
#         Storage = 0.
#         for i in 1:1000
#             Effective_Precipitation = rand(0.1:0.5:15)[1]
#             Potential_Evaporation = rand(0.1:0.2:5)[1]
#             Interception_Evaporation = rand(0:0.1:Potential_Evaporation * 0.5)[1]
#             Overlandflow, Preferentialflow, Soil_Evaporation, Soilstorage = soilstorage(Effective_Precipitation, Interception_Evaporation, Potential_Evaporation, Storage, beta, Ce, Ratio_Pref, Soilstoragecapacity)
#             Storage = Soilstorage
#             @test Storage <= Soilstoragecapacity
#         end
#     end
# end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
module Rules

import IterTools
import PlayingCards: Card, rankint

export cardvalue, validplays

function cardvalue(rankint::Integer)::Integer
    return min(rankint, 10)
end

function cardvalue(card::Card)::Integer
    return cardvalue(rankint(card))
end

function validhands(cards::Tuple{Vararg{Card}})
    return [ s for s in IterTools.subsets(cards, 4) ]
end

function validplays(prior::Tuple{Vararg{Card}}, cards::Tuple{Vararg{Card}})::Array{Card,1}
    if length(prior) == 0
        priortotal = 0
    else
        priortotal = sum(cardvalue(c) for c in prior)
    end
    return [ c for c in cards if cardvalue(c) ≤ (31-priortotal) ]
end

end  # module Rules