text:
using MemoryMutate

include("testutils.jl")

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

expr = :(a->b->c->d)

arr = zeros(Float64,10)
ptr = pointer(arr)
arr_ptr = convert(Ptr{Float64},pointer_from_objref(arr)) # is ≈64 Bytes below ptr
# because (?) in most architectures heap "grows upwards" where stack "grows downwards"
# Why stack grows down https://gist.github.com/cpq/8598782

ptr2 = ptr + (3-1)*sizeof(Float64)
view = unsafe_wrap(Array{Float64,1},ptr2,(1,); own=false)
view[1] = 2.0

function setat(ptr::Ptr{T},i::Int64,x::T) where T
    ptr2 = ptr + (i-1)*sizeof(T)
    view = unsafe_wrap(Array{T,1},ptr2,(1,); own=false)
    view[1] = x
end

arr = Float32[0f0,1f0,2f0,3f0,4f0]
arrptr = pointer(arr)
v = 9f0; @assert arr[2] != v; @mem arrptr[2] = v; @assert arr[2] == v;

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

        struct G2; x :: Float32;               end # struct G2 { const float x; };
        struct F2; x :: Float32; g ::     G2 ; end # struct F2 { const float x; const G2  g; };
        struct E2; x :: Float32; f ::     F2 ; end # struct E2 { const float x; const F2  g; };
        struct D2; x :: Float32; e :: Ptr{E2}; end # struct D2 {       float x;       E2  g; };
        struct C2; x :: Float32; d ::     D2 ; end # struct C2 { const float x; const D2& g; };
        struct B2; x :: Float32; c ::     C2 ; end # struct B2 { const float x; const C2& g; };
mutable struct A2; x :: Float32; b ::     B2 ; end # struct A2 {       float x;       B2& g; };

g2    =     G2(6.0f0)     # NOTE: we cannot change g2 and f2, since they aren't allocated on the heap
f2    =     F2(5.0f0,g2)  #
e2    = Ref(E2(4.0f0,f2)) # here, we copy f2 (and it's inner g2) to the heap -> these ones, e2[].f, e2[].f.g and e2[].f.g.x, can be changed
e2ptr = pointer_from_objref(e2)
d2    =     D2(3.0f0,e2ptr)
c2    =     C2(2.0f0,d2)
b2    =     B2(1.0f0,c2)
a2    =     A2(0.0f0,b2)

test201(a::A2,v::Float32) = @mem  a.b.c.d.e->f.g.x = v
test202(a::A2,v::G2)      = @mem  a.b.c.d.e->f.g = v
test203(a::A2,v::Float32) = @mem  a.b.c.d.e->f.x = v
test204(a::A2,v::F2)      = @mem  a.b.c.d.e->f = v
test205(a::A2,v::Float32) = @mem  a.b.c.d.e->x = v
test206(a::A2,v::Ptr{E2}) = @mem  a.b.c.d.e = v
test207(a::A2,v::Float32) = @mem  a.b.c.d.x = v
test208(a::A2,v::D2)      = @mem  a.b.c.d = v
test209(a::A2,v::Float32) = @mem  a.b.c.x = v
test210(a::A2,v::C2)      = @mem  a.b.c = v
test211(a::A2,v::Float32) = @mem  a.b.x = v
test212(a::A2,v::B2)      = @mem  a.b = v
test213(a::A2,v::Float32) = @mem  a.x = v
test214(a::A2,v::A2)      = @mem  a = v

v  = 0f0
v += 1f0;                                @assert unsafe_load(a2.b.c.d.e).f.g.x  != v;           test201(a2,v);           @assert unsafe_load(a2.b.c.d.e).f.g.x  == v;
v += 1f0;                                @assert unsafe_load(a2.b.c.d.e).f.g    != G2(v);       test202(a2,G2(v));       @assert unsafe_load(a2.b.c.d.e).f.g    == G2(v);
v += 1f0;                                @assert unsafe_load(a2.b.c.d.e).f.x    != v;           test203(a2,v);           @assert unsafe_load(a2.b.c.d.e).f.x    == v;
v += 1f0;                                @assert unsafe_load(a2.b.c.d.e).f      != F2(v,g2);    test204(a2,F2(v,g2));    @assert unsafe_load(a2.b.c.d.e).f      == F2(v,g2);
v += 1f0;                                @assert unsafe_load(a2.b.c.d.e).x      != v;           test205(a2,v);           @assert unsafe_load(a2.b.c.d.e).x      == v;
v += 1f0; e2b, e2bptr = mkptr(E2(v,f2)); @assert             a2.b.c.d.e         != e2bptr;      test206(a2,e2bptr);      @assert             a2.b.c.d.e         == e2bptr;
v += 1f0;                                @assert             a2.b.c.d.x         != v;           test207(a2,v);           @assert             a2.b.c.d.x         == v;
v += 1f0;                                @assert             a2.b.c.d           != D2(v,e2ptr); test208(a2,D2(v,e2ptr)); @assert             a2.b.c.d           == D2(v,e2ptr);
v += 1f0;                                @assert             a2.b.c.x           != v;           test209(a2,v);           @assert             a2.b.c.x           == v;
v += 1f0;                                @assert             a2.b.c             != C2(v,d2);    test210(a2,C2(v,d2));    @assert             a2.b.c             == C2(v,d2);
v += 1f0;                                @assert             a2.b.x             != v;           test211(a2,v);           @assert             a2.b.x             == v;
v += 1f0;                                @assert             a2.b               != B2(v,c2);    test212(a2,B2(v,c2));    @assert             a2.b               == B2(v,c2);
v += 1f0;                                @assert             a2.x               != v;           test213(a2,v);           @assert             a2.x               == v;
v += 1f0;                                @assert             a2                 != A2(v,b2);    test214(a2,A2(v,b2));  # @assert             a2                 == A2(v,b2); # # call by "sharing"

code_native(io,test201,(A2,Float32)); display_asm_stat_io(io) # (total = 3, movs = 2, mov = 1, vmov = 1)
code_native(io,test202,(A2,G2));      display_asm_stat_io(io) # (total = 6, movs = 4, mov = 3, vmov = 1)
code_native(io,test203,(A2,Float32)); display_asm_stat_io(io) # (total = 3, movs = 2, mov = 1, vmov = 1)
code_native(io,test204,(A2,F2));      display_asm_stat_io(io) # (total = 6, movs = 5, mov = 3, vmov = 2)
code_native(io,test205,(A2,Float32)); display_asm_stat_io(io) # (total = 4, movs = 2, mov = 1, vmov = 1)
code_native(io,test206,(A2,Ptr{E2})); display_asm_stat_io(io) # (total = 4, movs = 2, mov = 2, vmov = 0)
code_native(io,test207,(A2,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test208,(A2,D2));      display_asm_stat_io(io) # (total = 6, movs = 5, mov = 1, vmov = 4)
code_native(io,test209,(A2,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test210,(A2,C2));      display_asm_stat_io(io) # (total = 10, movs = 9, mov = 5, vmov = 4)
code_native(io,test211,(A2,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test212,(A2,B2));      display_asm_stat_io(io) # (total = 8, movs = 5, mov = 1, vmov = 4)
code_native(io,test213,(A2,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test214,(A2,A2));      display_asm_stat_io(io) # (total = 3, movs = 2, mov = 2, vmov = 0)

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

struct G3; x :: Float32;               end
struct F3; x :: Float32; g ::     G3 ; end
struct E3; x :: Float32; f ::     F3 ; end
struct D3; x :: Float32; e :: Ptr{E3}; end
struct C3; x :: Float32; d ::     D3 ; end
struct B3; x :: Float32; c ::     C3 ; end
struct A3; x :: Float32; b ::     B3 ; end

g3    =     G3(6.0f0)     # NOTE: we cannot change g2 and f2, since they aren't allocated on the heap
f3    =     F3(5.0f0,g3)  #
e3    = Ref(E3(4.0f0,f3)) # here, we copy f2 (and it's inner g2) to the heap -> these ones, e2[].f, e2[].f.g and e2[].f.g.x, can be changed
e3ptr = Ptr{E3}(pointer_from_objref(e3))
d3    =     D3(3.0f0,e3ptr)
c3    =     C3(2.0f0,d3)
b3    =     B3(1.0f0,c3)
a3    = Ref(A3(0.0f0,b3))
a3ptr = Ptr{A3}(pointer_from_objref(a3))

@assert unsafe_load(unsafe_load(a3ptr).b.c.d.e).f.g.x == @mem a3ptr->b.c.d.e->f.g.x
@assert unsafe_load(unsafe_load(a3ptr).b.c.d.e).f.g   == @mem a3ptr->b.c.d.e->f.g
@assert unsafe_load(unsafe_load(a3ptr).b.c.d.e).f     == @mem a3ptr->b.c.d.e->f
@assert unsafe_load(unsafe_load(a3ptr).b.c.d.e)       == @mem a3ptr->b.c.d.e[]
@assert             unsafe_load(a3ptr).b.c.d.e        == @mem a3ptr->b.c.d.e
@assert             unsafe_load(a3ptr).b.c.d          == @mem a3ptr->b.c.d
@assert             unsafe_load(a3ptr).b.c            == @mem a3ptr->b.c
@assert             unsafe_load(a3ptr).b              == @mem a3ptr->b
@assert             unsafe_load(a3ptr)                == @mem a3ptr[]

test301(a::Ptr{A3},v::Float32) = @mem  a->b.c.d.e->f.g.x = v
test302(a::Ptr{A3},v::G3)      = @mem  a->b.c.d.e->f.g = v
test303(a::Ptr{A3},v::Float32) = @mem  a->b.c.d.e->f.x = v
test304(a::Ptr{A3},v::F3)      = @mem  a->b.c.d.e->f = v
test305(a::Ptr{A3},v::Float32) = @mem  a->b.c.d.e->x = v
test306(a::Ptr{A3},v::Ptr{E3}) = @mem  a->b.c.d.e = v
test307(a::Ptr{A3},v::Float32) = @mem  a->b.c.d.x = v
test308(a::Ptr{A3},v::D3)      = @mem  a->b.c.d = v
test309(a::Ptr{A3},v::Float32) = @mem  a->b.c.x = v
test310(a::Ptr{A3},v::C3)      = @mem  a->b.c = v
test311(a::Ptr{A3},v::Float32) = @mem  a->b.x = v
test312(a::Ptr{A3},v::B3)      = @mem  a->b = v
test313(a::Ptr{A3},v::Float32) = @mem  a->x = v
test314(a::Ptr{A3},v::A3)      = @mem  a[] = v

v  = 0f0
v += 1f0;                                @assert unsafe_load(unsafe_load(a3ptr).b.c.d.e).f.g.x  != v;           test301(a3ptr,v);           @assert unsafe_load(unsafe_load(a3ptr).b.c.d.e).f.g.x  == v;
v += 1f0;                                @assert unsafe_load(unsafe_load(a3ptr).b.c.d.e).f.g    != G3(v);       test302(a3ptr,G3(v));       @assert unsafe_load(unsafe_load(a3ptr).b.c.d.e).f.g    == G3(v);
v += 1f0;                                @assert unsafe_load(unsafe_load(a3ptr).b.c.d.e).f.x    != v;           test303(a3ptr,v);           @assert unsafe_load(unsafe_load(a3ptr).b.c.d.e).f.x    == v;
v += 1f0;                                @assert unsafe_load(unsafe_load(a3ptr).b.c.d.e).f      != F3(v,g3);    test304(a3ptr,F3(v,g3));    @assert unsafe_load(unsafe_load(a3ptr).b.c.d.e).f      == F3(v,g3);
v += 1f0;                                @assert unsafe_load(unsafe_load(a3ptr).b.c.d.e).x      != v;           test305(a3ptr,v);           @assert unsafe_load(unsafe_load(a3ptr).b.c.d.e).x      == v;
v += 1f0; e3b, e3bptr = mkptr(E3(v,f3)); @assert             unsafe_load(a3ptr).b.c.d.e         != e3bptr;      test306(a3ptr,e3bptr);      @assert             unsafe_load(a3ptr).b.c.d.e         == e3bptr;
v += 1f0;                                @assert             unsafe_load(a3ptr).b.c.d.x         != v;           test307(a3ptr,v);           @assert             unsafe_load(a3ptr).b.c.d.x         == v;
v += 1f0;                                @assert             unsafe_load(a3ptr).b.c.d           != D3(v,e3ptr); test308(a3ptr,D3(v,e3ptr)); @assert             unsafe_load(a3ptr).b.c.d           == D3(v,e3ptr);
v += 1f0;                                @assert             unsafe_load(a3ptr).b.c.x           != v;           test309(a3ptr,v);           @assert             unsafe_load(a3ptr).b.c.x           == v;
v += 1f0;                                @assert             unsafe_load(a3ptr).b.c             != C3(v,d3);    test310(a3ptr,C3(v,d3));    @assert             unsafe_load(a3ptr).b.c             == C3(v,d3);
v += 1f0;                                @assert             unsafe_load(a3ptr).b.x             != v;           test311(a3ptr,v);           @assert             unsafe_load(a3ptr).b.x             == v;
v += 1f0;                                @assert             unsafe_load(a3ptr).b               != B3(v,c3);    test312(a3ptr,B3(v,c3));    @assert             unsafe_load(a3ptr).b               == B3(v,c3);
v += 1f0;                                @assert             unsafe_load(a3ptr).x               != v;           test313(a3ptr,v);           @assert             unsafe_load(a3ptr).x               == v;
v += 1f0; a3b, a3bptr = mkptr(A3(v,b3)); @assert             unsafe_load(a3ptr)                 != A3(v,b3);    test314(a3ptr,A3(v,b3));    @assert             unsafe_load(a3ptr)                 == A3(v,b3);

code_native(io,test301,(Ptr{A3},Float32)); display_asm_stat_io(io) # (total = 3, movs = 2, mov = 1, vmov = 1)
code_native(io,test302,(Ptr{A3},G3));      display_asm_stat_io(io) # (total = 6, movs = 4, mov = 3, vmov = 1)
code_native(io,test303,(Ptr{A3},Float32)); display_asm_stat_io(io) # (total = 3, movs = 2, mov = 1, vmov = 1)
code_native(io,test304,(Ptr{A3},F3));      display_asm_stat_io(io) # (total = 6, movs = 5, mov = 3, vmov = 2)
code_native(io,test305,(Ptr{A3},Float32)); display_asm_stat_io(io) # (total = 4, movs = 2, mov = 1, vmov = 1)
code_native(io,test306,(Ptr{A3},Ptr{E3})); display_asm_stat_io(io) # (total = 4, movs = 2, mov = 2, vmov = 0)
code_native(io,test307,(Ptr{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test308,(Ptr{A3},D3));      display_asm_stat_io(io) # (total = 6, movs = 5, mov = 1, vmov = 4)
code_native(io,test309,(Ptr{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test310,(Ptr{A3},C3));      display_asm_stat_io(io) # (total = 10, movs = 9, mov = 5, vmov = 4)
code_native(io,test311,(Ptr{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test312,(Ptr{A3},B3));      display_asm_stat_io(io) # (total = 8, movs = 5, mov = 1, vmov = 4)
code_native(io,test313,(Ptr{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test314,(Ptr{A3},A3));      display_asm_stat_io(io) # (total = 7, movs = 5, mov = 3, vmov = 2) # TODO: got (total = 11, movs = 9, mov = 5, vmov = 4)


text:
module vSmartMOM

using LinearAlgebra                # For linear algebra routines
using NCDatasets                   # For reading netcdf atmospheric profiles
using ProgressMeter                # Showing progress in for loops
using Distributions                # Distributions of aerosols
using DelimitedFiles               # For reading ASCII files
using ..Scattering                 # Use scattering module
using ..Absorption                 # Use absorption module
using ...RadiativeTransfer         # Use parent RadiativeTransfer module
using ...Architectures             # Use Architectures module

using CUDA                         # GPU CuArrays and functions
using KernelAbstractions           # Abstracting code for CPU/GPU
using KernelAbstractions.Extras

using FastGaussQuadrature          # Computes quadrature points (Gauss, legendre, Radau,...)
using TimerOutputs                 # For timing sections of the code
using StatsBase                    # Fit statistics for truncation
using Parameters                   # For keyword arguments in structs
using DocStringExtensions          # For documenting
using YAML                         # For reading properties files
using ForwardDiff                  # Automatic Differentiation
using NNlib                        # For batched multiplications
import NNlib.batched_mul           # Required to overwrite batched_mul for Duals

# More threads in LA wasn't really helpful, can be turned off here:
# LinearAlgebra.BLAS.set_num_threads(1)

# Constants and Types
include("Constants/constants.jl")         # Scientific constants
include("Types/types.jl")                 # All custom types for this module

# Solvers
include("Solvers/elemental.jl")           # Elemental
include("Solvers/doubling.jl")            # Doubling
include("Solvers/interaction.jl")         # Interaction
include("Solvers/rt_kernel.jl")           # Handle Core RT (Elemental/Doubling/Interaction)
include("Solvers/postprocessing_vza.jl")  # Postprocess (Azimuthal Weighting)
include("Solvers/rt_run.jl")              # Starting point for RT

# GPU
include("GPU/gpu_batched.jl")             # Batched operations
include("GPU/CUDA_getri.jl")              # Custom getri_strided_batched!

# Utilities / Helper Functions
include("Utils/atmo_prof.jl")             # Helper Functions for Hanling Atmospheric Profiles
include("Utils/rt_utils.jl")              # Miscellaneous Utility Functions
include("Utils/rt_streams.jl")            # Set streams before RT
include("Utils/forwardDiff_tools.jl")     # Helpers for Forward Differentiation
include("Utils/model_parameters.jl")      # Handling Model Parameters

# Surfaces
include("Surface/lambertian_surface.jl")  # Lambertian Surface

# Functions to export
export parameters_from_yaml,              # Getting parameters from a file
       model_from_parameters,             # Converting the parameters to model
       rt_run,                            # Run the RT code
       default_parameters                 # Set of default parameters

end

ext: jl


text:
# Script for generating upstream MHD fluctuations.
#
# Hongyang Zhou, hyzhou@umich.edu

using PyPlot, SpaceAnalysis

## Background values

n0 = 2e6 # []
T0 = 5e5
V0 = [-600., 0.0, 0.0]*1e3 # [m/s]
B0 = [5., 0., 0.]*1e-9 # [T]

varBG = BackgroundVariable(n0, T0, V0..., B0...)

fsignal = 1/150 # perturbation frequency, [Hz]
signal  = :fast # perturbation type
dv      = 5e4   # perturbation magnitude in velocity, [m/s]
dn      = 1     # perturbation magnitude in density, [amu/cc]
dir     = "y"   # perturbed direction(s)
fsample = 2     # sampling frequency, [Hz]
tstart  = 300   # [s]
tend    = 750   # [s]

s = generate_signal(varBG; fsignal, signal, dv, fsample, tstart, tend, dir)

save_signal(s)

## Visualization

fig, ax1 = plt.subplots(nrows=1, ncols=1)

if signal == :alfven
   color = "tab:red"
   ax1.set_xlabel("time (s)")
   ax1.set_ylabel("By (nT)"; color)
   ax1.plot(s.t, s.By*1e9; color)
   ax1.tick_params(axis="y", labelcolor=color)

   ax2 = ax1.twinx()
   color = "tab:blue"
   ax2.set_ylabel("Uy (km/s)"; color)
   ax2.plot(s.t, s.Vy/1e3, "--"; color)
   ax2.tick_params(axis="y", labelcolor=color)
elseif signal == :fast
   color = "tab:red"
   ax1.set_xlabel("time (s)")
   ax1.set_ylabel("density (amu/cc)"; color)
   ax1.plot(s.t, s.n*1e-6; color)
   ax1.tick_params(axis="y", labelcolor=color)
end


text:
"""
senorge_info()

Get coordinate information about senorge grids
"""
function senorge_info()

nrows, ncols = 1550, 1195

xllcorner, yllcorner = -75000, 6450000

senorge_ind = reshape(collect(1:nrows*ncols), ncols, nrows)'

xvec = collect(xllcorner:1000:(xllcorner+(ncols-1)*1000))'

xcoord = repmat(xvec, nrows, 1)

yvec = collect(yllcorner+(nrows-1)*1000:-1000:yllcorner)

ycoord = repmat(yvec, 1, ncols)

return senorge_ind, xcoord, ycoord

end

"""
read_esri_raster(file_name)

Read esri raster file.
"""
function read_esri_raster(file_name)

lines = readlines(file_name)

res = Dict()

for irow in 1:6
    str_splitted = split(lines[irow])
    res[str_splitted[1]] = parse(Int64, str_splitted[2])
end

raster = readdlm(file_name, ' ', skipstart=6)

res["data"] = convert(Array{Float64,2}, raster[:, 1:end-1])

return res

end

"""
raster2vec(raster)

Convert a esri raster to a vector and replace missing data with not-a-number.
"""
function raster2vec(raster)

data = raster["data"]
data[data .== raster["NODATA_value"]] .= NaN
data = permutedims(data)
vec = data[:]

return vec

end

"""
read_bil(file_name)

Read binary file with senorge data.
"""
function read_bil(file_name)

data = zeros(Float64, 1550*1195)

read_bil!(data, file_name)

return data

end

function read_bil!(data, file_name)

@assert length(data) == 1550*1195

f = open(file_name, "r")
i = 1
while ~eof(f)
    data[i] = read(f, Int16)
    i += 1
end
close(f)

end

max_forks_repo_name: jmgnve/NveData


text:
using StochasticRounding
using Test

include("bfloat16sr.jl")
include("float16sr.jl")
include("float32sr.jl")
include("seeding.jl")


text:
# Li, Lin, "Accelerated Proximal Gradient Methods for Nonconvex Programming",
# Proceedings of NIPS 2015 (2015).

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

"""
    LiLinIteration(; <keyword-arguments>)

Instantiate the nonconvex accelerated proximal gradient method by Li and Lin
(see Algorithm 2 in [1]) for solving optimization problems of the form

    minimize f(x) + g(x),

where `f` is smooth.

# Arguments
- `x0`: initial point.
- `f=Zero()`: smooth objective term.
- `g=Zero()`: proximable objective term.
- `Lf=nothing`: Lipschitz constant of the gradient of f.
- `gamma=nothing`: stepsize to use, defaults to `1/Lf` if not set (but `Lf` is).

# References
- [1] Li, Lin, "Accelerated Proximal Gradient Methods for Nonconvex Programming",
Proceedings of NIPS 2015 (2015).
"""

Base.@kwdef struct LiLinIteration{R,C<:Union{R,Complex{R}},Tx<:AbstractArray{C},Tf,Tg}
    f::Tf = Zero()
    g::Tg = Zero()
    x0::Tx
    Lf::Maybe{R} = nothing
    gamma::Maybe{R} = Lf === nothing ? nothing : (1 / Lf)
    adaptive::Bool = false
    delta::R = real(eltype(x0))(1e-3)
    eta::R = real(eltype(x0))(0.8)
end

Base.IteratorSize(::Type{<:LiLinIteration}) = Base.IsInfinite()

mutable struct LiLinState{R<:Real,Tx}
    x::Tx             # iterate
    y::Tx             # extrapolated point
    f_y::R            # value of f at y
    grad_f_y::Tx      # gradient of f at y
    # TODO: *two* gammas should be used in general, one for y and one for x
    gamma::R          # stepsize parameter of forward and backward steps
    y_forward::Tx     # forward point at y
    z::Tx             # forward-backward point
    g_z::R            # value of nonsmooth term at z
    res::Tx           # fixed-point-residual (at y)
    theta::R          # auxiliary sequence to compute extrapolated points
    F_average::R      # moving average of objective values
    q::R              # auxiliary sequence to compute moving average
end

function Base.iterate(iter::LiLinIteration{R}) where {R}
    y = copy(iter.x0)
    grad_f_y, f_y = gradient(iter.f, y)

    # TODO: initialize gamma if not provided
    # TODO: authors suggest Barzilai-Borwein rule?
    # TODO: *two* gammas should be used in general, one for y and one for x

    # compute initial forward-backward step
    y_forward = y - iter.gamma .* grad_f_y
    z, g_z = prox(iter.g, y_forward, iter.gamma)

    Fy = f_y + iter.g(y)

    @assert isfinite(Fy) "initial point must be feasible"

    # compute initial fixed-point residual
    res = y - z

    state = LiLinState(
        copy(iter.x0), y, f_y, grad_f_y, iter.gamma,
        y_forward, z, g_z, res, R(1), Fy, R(1),
    )

    return state, state
end

function Base.iterate(
    iter::LiLinIteration{R},
    state::LiLinState{R,Tx},
) where {R,Tx}
    # TODO: backtrack gamma at y

    Fz = iter.f(state.z) + state.g_z

    theta1 = (R(1) + sqrt(R(1) + 4 * state.theta^2)) / R(2)

    if Fz <= state.F_average - iter.delta * norm(state.res)^2
        case = 1
    else
        # TODO: re-use available space in state?
        # TODO: backtrack gamma at x
        grad_f_x, f_x = gradient(iter.f, x)
        x_forward = state.x - state.gamma .* grad_f_x
        v, g_v = prox(iter.g, x_forward, state.gamma)
        Fv = iter.f(v) + g_v
        case = Fz <= Fv ? 1 : 2
    end

    if case == 1
        state.y .= state.z .+ ((state.theta - R(1)) / theta1) .* (state.z .- state.x)
        state.x, state.z = state.z, state.x
        Fx = Fz
    elseif case == 2
        state.y .=
            state.z .+ (state.theta / theta1) .* (state.z .- v) .+
            ((state.theta - R(1)) / theta1) .* (v .- state.x)
        state.x = v
        Fx = Fv
    end

    state.f_y = gradient!(state.grad_f_y, iter.f, state.y)
    state.y_forward .= state.y .- state.gamma .* state.grad_f_y
    state.g_z = prox!(state.z, iter.g, state.y_forward, state.gamma)

    state.res .= state.y - state.z

    state.theta = theta1

    # NOTE: the following can be simplified
    q1 = iter.eta * state.q + 1
    state.F_average = (iter.eta * state.q * state.F_average + Fx) / q1
    state.q = q1

    return state, state
end

# Solver

struct LiLin{R, K}
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int
    kwargs::K
end

function (solver::LiLin)(x0; kwargs...)
    stop(state::LiLinState) = norm(state.res, Inf) / state.gamma <= solver.tol
    disp((it, state)) =
        @printf("%5d | %.3e | %.3e\n", it, state.gamma, norm(state.res, Inf) / state.gamma)
    iter = LiLinIteration(; x0=x0, solver.kwargs..., kwargs...)
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose
        iter = tee(sample(iter, solver.freq), disp)
    end
    num_iters, state_final = loop(iter)
    return state_final.z, num_iters
end

LiLin(; maxit=10_000, tol=1e-8, verbose=false, freq=100, kwargs...) =
    LiLin(maxit, tol, verbose, freq, kwargs)


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".

function _co_refract(old_alt::T, altitude::T, pressure::T, temperature::T,
                    epsilon::T, to_observe::Bool) where {T<:AbstractFloat}

    if isnan(temperature)
        if altitude > 11000
            temperature = T(211.5)
        else
            temperature = 283 - 0.0065 * altitude
        end
    end

    if isnan(pressure)
        pressure = 1010 * (1 - 6.5/288000 * altitude) ^ 5.255
    end

    if !to_observe
        aout = old_alt - co_refract_forward(old_alt, pressure, temperature)
    else
        cur = old_alt + co_refract_forward(old_alt, pressure, temperature)
        while true
            last = cur
            cur = old_alt + co_refract_forward(last, pressure, temperature)
            if abs(last - cur) * 3600 < epsilon
                break
            end
        end
        aout = cur
    end
    return aout
end
"""
    co_refract(old_alt[, altitude=0, pressure=NaN, temperature=NaN,
               epsilon=0.25, to_observe=false]) -> aout

### Purpose ###

Calculate correction to altitude due to atmospheric refraction.

### Explanation ##

Because the index of refraction of air is not precisely 1.0, the atmosphere
bends all incoming light, making a star or other celestial object appear at
a slightly different altitude (or elevation) than it really is.  It is
important to understand the following definitions:

* Observed Altitude: The altitude that a star is seen to be, with a telescope.
  This is where it appears in the sky. This is should be always greater than the
  apparent altitude.

* Apparent Altitude: The altitude that a star would be at, if ~there were no
  atmosphere~ (sometimes called the "true" altitude). This is usually calculated
  from an object's celestial coordinates. Apparent altitude should always be
  smaller than the observed altitude.

Thus, for example, the Sun's apparent altitude when you see it right on the
horizon is actually -34 arcminutes.

This program uses a couple of simple formulae to estimate the effect for most
optical and radio wavelengths. Typically, you know your observed altitude
(from an observation), and want the apparent altitude. To go the other way,
this program uses an iterative approach.

### Arguments ###

* `old_alt`: observed altitude in degrees. If `to_observe` is set to true,
  this should be apparent altitude
* `altitude` (optional): the height of the observing location, in meters. This is
  only used to determine an approximate temperature and pressure, if these are not
  specified separately. Default is 0 i.e. sea level
* `pressure` (optional): the pressure at the observing location, in millibars.
  Default is NaN
* `temperature` (optional): the temperature at the observing location, in Kelvins.
  Default is NaN
* `epsilon` (optional): the accuracy to obtain, in arcseconds. If `to_observe` is
  true, then it will be calculated. Default is 0.25 arcseconds
* `to_observe` (optional boolean keyword): if set to true, it is assumed that `old_alt`
  has apparent altitude as its input and the observed altitude will be found

### Output ###

* `aout`: apparent altitude, in degrees. Observed altitude is returned if `to_observe`
  is set to true

### Example ###

The lower limb of the Sun is observed to have altitude of 0d 30'. Calculate the the true
(i.e. apparent) altitude of the Sun's lower limb using mean  conditions of air pressure and
temperature.

```jldoctest
julia> using AstroLib

julia> co_refract(0.5)
0.02584736873098442
```

### Notes ###

If altitude is set but the temperature or pressure is not, the program will make an
intelligent guess for the temperature and pressure.

#### Wavelength Dependence ####

This correction is 0 at zenith, about 1 arcminute at 45 degrees, and 34 arcminutes at the
horizon for optical wavelengths. The correction is non-negligible at all wavelengths,
but is not very easily calculable.
These formulae assume a wavelength of 550 nm, and will be accurate to about 4 arcseconds
for all visible wavelengths, for elevations of 10 degrees and higher.
Amazingly, they are also accurate for radio frequencies less than ~ 100 GHz.

#### References ####

* Meeus, Astronomical Algorithms, Chapter 15.
* Explanatory Supplement to the Astronomical Almanac, 1992.
* Methods of Experimental Physics, Vol 12 Part B, Astrophysics, Radio Telescopes, Chapter 2.5,
  "Refraction Effects in the Neutral Atmosphere", by R.K. Crane.

Code of this function is based on IDL Astronomy User's Library.
"""
co_refract(old_alt::Real, altitude::Real=0, pressure::Real=NaN, temperature::Real=NaN,
           epsilon::Real=0.25; to_observe::Bool=false) =
               _co_refract(promote(float(old_alt), float(altitude), float(pressure),
                          float(temperature), float(epsilon))..., to_observe)

"""
    co_refract_forward(alt, pre, temp) -> ref

### Purpose ###

A function used by [`co_refract`](@ref) to find apparent (or observed) altitude

### Arguments ###

* `alt`: the observed (or apparent) altitude, in degrees
* `pre`: pressure, in millibars
* `temp`: temperature, in Kelvins

### Output ###

* `ref`: the atmospheric refraction, in minutes of arc

### Notes ###

The atmospheric refraction is calculated by Saemundsson's formula

Code of this function is based on IDL Astronomy User's Library.
"""
function co_refract_forward(alt::Real, pre::Real, temp::Real)

    if alt < 15
        ref =3.569*@evalpoly(alt, 0.1594, 0.0196, 0.00002)/@evalpoly(alt, 1, 0.505, 0.0845)
    else
        ref = 0.0166667 / tand((alt + 7.31 / (alt + 4.4)))
    end
    ref *= pre * 283 / (1010 * temp)
    return ref
end


text:
using PyPlot, BinningAnalysis

function performance_test(N = 28, start=1)
    timings = Float64[]

    for i in start:N
        # NOTE using the default size here actually performs better
        BA = LogBinner(i+1)
        xs = rand(2^i)
        GC.gc()

        t0 = time_ns()
        for x in xs
            push!(BA, x)
        end
        t = time_ns() - t0

        println("2^$i done in $(t)ns.")
        push!(timings, t / length(xs))
    end

    start:N, timings
end

performance_test(1)
xs, ys = performance_test(26)

begin
    figure()
    plot(xs, ys, "r+:")
    xlabel("log(Values pushed)")
    ylabel("Average time per push [ns]")
    ylim(0, 100)
end

# Memory: 33B * N where N ≥ ⌈log(K+1)⌉ with K values pushed
# Time:   10-40ns per push, depending on luck and how I count


text:
# This file includes callbacks for dynamical networks.

"""
    $TYPEDSIGNATURES

Returns a `PCMCallback`. `nbits` is the number of bits, `ti` is the initial time, `tbit` is the bit time, `high` is the high level of coupling and `low` is the low level of coupling.
"""
function PCMCallback(nbits::Int = 10, tbit::Real = 1.0;
    ti::Real = 0.0, high::Real = 1.0, low::Real = 0.0001, src::Int = 1, dst::Int = 2)
    tf = tbit * nbits
    tstops = ti:tbit:tf-tbit
    condition(u, t, integrator) = integrator.t ∈ tstops
    function action!(integrator)
        net, bits = integrator.sol.prob.p
        popfirst!(bits) ? add_edge!(net.graph, src, dst, high) : add_edge!(net.graph, src, dst, low)
        # integrator.u += rand(length(integrator.u))  # Perturbation to desynchronize
    end
    DiscreteCallback(condition, action!)
end


text:
function git_make_commit(; commit_message::String)
    result = try
        cmd = `git commit -m $(commit_message)`
        p = pipeline(cmd; stdout=stdout, stderr=stderr)
        success(p)
    catch
        false
    end
    return result
end


text:
print("Loading module\n")

using OpenCV
using Test

test_dir = joinpath(ENV["OPENCV_TEST_DATA_PATH"], "cv")

include("test_feature2d.jl")
include("test_imgproc.jl")
include("test_objdetect.jl")

exit(0)


text:
include("TwoPartCost.jl")
include("ThreePartCost.jl")
include("MultiStartCost.jl")
include("Area.jl")
include("LoadZone.jl")
include("Bus.jl")
include("Arc.jl")
include("Line.jl")
include("MonitoredLine.jl")
include("PhaseShiftingTransformer.jl")
include("TapTransformer.jl")
include("Transformer2W.jl")
include("HVDCLine.jl")
include("VSCDCLine.jl")
include("InterruptibleLoad.jl")
include("FixedAdmittance.jl")
include("PowerLoad.jl")
include("HydroEnergyReservoir.jl")
include("HydroDispatch.jl")
include("RenewableDispatch.jl")
include("RenewableFix.jl")
include("ThermalStandard.jl")
include("ThermalMultiStart.jl")
include("GenericBattery.jl")
include("StaticReserve.jl")
include("ReserveDemandCurve.jl")
include("VariableReserve.jl")
include("AGC.jl")
include("Transfer.jl")
include("AVRFixed.jl")
include("AVRSimple.jl")
include("AVRTypeI.jl")
include("AVRTypeII.jl")
include("BaseMachine.jl")
include("RoundRotorMachine.jl")
include("SalientPoleMachine.jl")
include("AndersonFouadMachine.jl")
include("FullMachine.jl")
include("MarconatoMachine.jl")
include("OneDOneQMachine.jl")
include("SimpleAFMachine.jl")
include("SimpleFullMachine.jl")
include("SimpleMarconatoMachine.jl")
include("PSSFixed.jl")
include("PSSSimple.jl")
include("SingleMass.jl")
include("FiveMassShaft.jl")
include("TGFixed.jl")
include("TGTypeI.jl")
include("TGTypeII.jl")
include("AverageConverter.jl")
include("FixedDCSource.jl")
include("LCLFilter.jl")
include("LCFilter.jl")
include("KauraPLL.jl")
include("VirtualInertia.jl")
include("ReactivePowerDroop.jl")
include("CurrentControl.jl")
include("Source.jl")

export get_Ae
export get_Be
export get_D
export get_D_12
export get_D_23
export get_D_34
export get_D_45
export get_D_ex
export get_D_hp
export get_D_ip
export get_D_lp
export get_H
export get_H_ex
export get_H_hp
export get_H_ip
export get_H_lp
export get_K0
export get_K_d
export get_K_ex
export get_K_hp
export get_K_i
export get_K_ip
export get_K_lp
export get_K_p
export get_K_ω
export get_Ka
export get_Ke
export get_Kf
export get_Kv
export get_L_1d
export get_L_1q
export get_L_ad
export get_L_aq
export get_L_d
export get_L_f1d
export get_L_ff
export get_L_q
export get_P_max
export get_P_min
export get_P_ref
export get_R
export get_R_1d
export get_R_1q
export get_R_f
export get_Se
export get_T1
export get_T2
export get_T3
export get_T4
export get_T5
export get_T_AA
export get_Ta
export get_Tc
export get_Td0_p
export get_Td0_pp
export get_Te
export get_Tf
export get_Tq0_p
export get_Tq0_pp
export get_Tr
export get_Ts
export get_V_pss
export get_V_ref
export get_Vf
export get_Vr_max
export get_Vr_min
export get_X_th
export get_Xd
export get_Xd_p
export get_Xd_pp
export get_Xl
export get_Xq
export get_Xq_p
export get_Xq_pp
export get_Y
export get_activepower
export get_activepower_flow
export get_activepowerlimits
export get_activepowerlimits_from
export get_activepowerlimits_to
export get_angle
export get_anglelimits
export get_arc
export get_area
export get_available
export get_b
export get_basepower
export get_basevoltage
export get_bias
export get_bus
export get_bustype
export get_capacity
export get_cf
export get_delta_t
export get_dynamic_injector
export get_efficiency
export get_energy
export get_eq_p
export get_ext
export get_fixed
export get_flowlimits
export get_from
export get_fuel
export get_inflow
export get_initial_ace
export get_initial_storage
export get_inputactivepowerlimits
export get_internal
export get_inv_d_fluxlink
export get_inv_q_fluxlink
export get_inverter_firing_angle
export get_inverter_taplimits
export get_inverter_xrc
export get_kad
export get_kd
export get_kffi
export get_kffv
export get_ki_pll
export get_kic
export get_kiv
export get_kp_pll
export get_kpc
export get_kpv
export get_kq
export get_kω
export get_lf
export get_lg
export get_load_response
export get_load_zone
export get_loss
export get_lv
export get_maxactivepower
export get_maxreactivepower
export get_model
export get_must_run
export get_n_states
export get_no_load
export get_number
export get_op_cost
export get_outputactivepowerlimits
export get_power_trajectory
export get_powerfactor
export get_primaryshunt
export get_primemover
export get_r
export get_ramplimits
export get_rate
export get_rating
export get_reactivepower
export get_reactivepower_flow
export get_reactivepowerlimits
export get_reactivepowerlimits_from
export get_reactivepowerlimits_to
export get_rectifier_firing_angle
export get_rectifier_taplimits
export get_rectifier_xrc
export get_requirement
export get_rf
export get_rg
export get_rv
export get_s_rated
export get_services
export get_shutdn
export get_start_time_limits
export get_start_types
export get_startup
export get_states
export get_status
export get_storage_capacity
export get_tap
export get_time_at_status
export get_timeframe
export get_timelimits
export get_to
export get_v_rated
export get_variable
export get_voltage
export get_voltagelimits
export get_x
export get_α
export get_γd
export get_γq
export get_τ_max
export get_τ_min
export get_ω_lp
export get_ωad
export get_ωf
export set_Ae!
export set_Be!
export set_D!
export set_D_12!
export set_D_23!
export set_D_34!
export set_D_45!
export set_D_ex!
export set_D_hp!
export set_D_ip!
export set_D_lp!
export set_H!
export set_H_ex!
export set_H_hp!
export set_H_ip!
export set_H_lp!
export set_K0!
export set_K_d!
export set_K_ex!
export set_K_hp!
export set_K_i!
export set_K_ip!
export set_K_lp!
export set_K_p!
export set_K_ω!
export set_Ka!
export set_Ke!
export set_Kf!
export set_Kv!
export set_L_1d!
export set_L_1q!
export set_L_ad!
export set_L_aq!
export set_L_d!
export set_L_f1d!
export set_L_ff!
export set_L_q!
export set_P_max!
export set_P_min!
export set_P_ref!
export set_R!
export set_R_1d!
export set_R_1q!
export set_R_f!
export set_Se!
export set_T1!
export set_T2!
export set_T3!
export set_T4!
export set_T5!
export set_T_AA!
export set_Ta!
export set_Tc!
export set_Td0_p!
export set_Td0_pp!
export set_Te!
export set_Tf!
export set_Tq0_p!
export set_Tq0_pp!
export set_Tr!
export set_Ts!
export set_V_pss!
export set_V_ref!
export set_Vf!
export set_Vr_max!
export set_Vr_min!
export set_X_th!
export set_Xd!
export set_Xd_p!
export set_Xd_pp!
export set_Xl!
export set_Xq!
export set_Xq_p!
export set_Xq_pp!
export set_Y!
export set_activepower!
export set_activepower_flow!
export set_activepowerlimits!
export set_activepowerlimits_from!
export set_activepowerlimits_to!
export set_angle!
export set_anglelimits!
export set_arc!
export set_area!
export set_available!
export set_b!
export set_basepower!
export set_basevoltage!
export set_bias!
export set_bus!
export set_bustype!
export set_capacity!
export set_cf!
export set_delta_t!
export set_dynamic_injector!
export set_efficiency!
export set_energy!
export set_eq_p!
export set_ext!
export set_fixed!
export set_flowlimits!
export set_from!
export set_fuel!
export set_inflow!
export set_initial_ace!
export set_initial_storage!
export set_inputactivepowerlimits!
export set_internal!
export set_inv_d_fluxlink!
export set_inv_q_fluxlink!
export set_inverter_firing_angle!
export set_inverter_taplimits!
export set_inverter_xrc!
export set_kad!
export set_kd!
export set_kffi!
export set_kffv!
export set_ki_pll!
export set_kic!
export set_kiv!
export set_kp_pll!
export set_kpc!
export set_kpv!
export set_kq!
export set_kω!
export set_lf!
export set_lg!
export set_load_response!
export set_load_zone!
export set_loss!
export set_lv!
export set_maxactivepower!
export set_maxreactivepower!
export set_model!
export set_must_run!
export set_n_states!
export set_no_load!
export set_number!
export set_op_cost!
export set_outputactivepowerlimits!
export set_power_trajectory!
export set_powerfactor!
export set_primaryshunt!
export set_primemover!
export set_r!
export set_ramplimits!
export set_rate!
export set_rating!
export set_reactivepower!
export set_reactivepower_flow!
export set_reactivepowerlimits!
export set_reactivepowerlimits_from!
export set_reactivepowerlimits_to!
export set_rectifier_firing_angle!
export set_rectifier_taplimits!
export set_rectifier_xrc!
export set_requirement!
export set_rf!
export set_rg!
export set_rv!
export set_s_rated!
export set_services!
export set_shutdn!
export set_start_time_limits!
export set_start_types!
export set_startup!
export set_states!
export set_status!
export set_storage_capacity!
export set_tap!
export set_time_at_status!
export set_timeframe!
export set_timelimits!
export set_to!
export set_v_rated!
export set_variable!
export set_voltage!
export set_voltagelimits!
export set_x!
export set_α!
export set_γd!
export set_γq!
export set_τ_max!
export set_τ_min!
export set_ω_lp!
export set_ωad!
export set_ωf!


text:
abstract type Weights{F<:Real} <: AbstractVector{F} end

@inline (u::Weights{F})(i::Integer) where {F} = u[i]
@inline (u::Weights{F})(i::Integer, b::Bool) where {F} = u[i, b]

"""
Uniform one weights.

## Example
```jldoctest
julia> w = GraphIdx.Ones{Int}(); w[3]
1

julia> w[5] = -2; w[5]
1
```
"""
struct Ones{F<:Real} <: Weights{F} end

@inline Base.getindex(::Ones{F}, ::Integer) where {F} = one(F)
@inline Base.getindex(::Ones{F}, ::Integer, ::Bool) where {F} = one(F)
@inline Base.setindex!(::Ones{F}, ::F, ::Integer) where {F} = one(F)
Base.collect(w::Ones{F}, n::Integer) where {F} = ones(F, n)

"""
Represent weighting that is indexable and callable where each index has the same value.

Usually, weights are assumed to be positive!

## Example
```jldoctest
julia> w = GraphIdx.Const(5.4);

julia> w[3]
5.4

julia> w[5] = 2.1; w[5]
5.4
```
"""
struct Const{F<:Real} <: Weights{F}
    w::F
end

@inline Base.getindex(c::Const{F}, ::Integer) where {F} = c.w
@inline Base.getindex(c::Const{F}, ::Integer, ::Bool) where {F} = c.w
@inline Base.setindex!(c::Const{F}, ::F, ::Integer) where {F} = c.w
Base.collect(w::Const{F}, n::Integer) where {F} = fill(w.w, n)

"""
Weighting of nodes where every weight can be different.
The actual value of a node can be accessed via call - or index syntax.

## Example
```jldoctest
julia> w = GraphIdx.Vec([0.1, 0.2, 0.3]);

julia> w[3]
0.3

julia> w[1]
0.1

julia> w[1] = 0; w[1]
0.0
```
"""
struct Vec{F<:Real} <: Weights{F}
    a::Vector{F}
end

@inline Base.getindex(v::Vec{F}, i::Integer) where {F} = v.a[i]
@inline Base.getindex(v::Vec{F}, i::Integer, ::Bool) where {F} = v.a[i]
@inline Base.setindex!(v::Vec{F}, c::W, i::Integer) where {F, W} = (v.a[i] = F(c))
@inline Base.collect(w::Vec{F}, n::Integer) where {F} =
    (@assert n == length(w.a) "n=$n != $(length(w.a))"; w.a)

"""
Signed vector (e.g. for directed edge weights)
"""
struct SVec{F<:Real} <: Weights{F}
    a::Array{F, 2}
end

SVec(v1::Vector{F}, v2::Vector{F}) where {F} = SVec(cat(v1, v2, dims=2))
SVec(v::NTuple{2, Vector{F}}) where {F} = SVec(cat(v..., dims=2))

@inline Base.getindex(v::SVec{F}, i::Integer) where {F} = v.a[i, 1]
@inline Base.getindex(v::SVec{F}, i::Integer, b::Bool) where {F} = v.a[i, Int(b) + 1]
@inline Base.setindex!(v::SVec{F}, c::W, i::Integer, b::Bool) where {F, W} =
    (v.a[i, Int(b) + 1] = F(c))
@inline function Base.collect(w::SVec{F}, n::Integer) where {F}
    @assert n == size(w.a, 1) "n=$n != $(size(w.a, 1))"
    @assert 2 == size(w.a, 2) "$(size(w.a, 2))"
    w.a
end

Weights(::Type{F}) where {F<:Real} = Ones{F}()
Weights(w::F) where {F<:Real} = w == F(1.0) ? Ones{F}() : Const(w)
Weights(w::Vector{F}) where {F<:Real} = Vec(w)
Weights(w::Array{F, 2}) where {F<:Real} = SVec(w)
Weights(w::NTuple{2, Vector{F}}) where {F<:Real} = SVec(cat(w..., dims=2))

Base.similar(::Ones{F}, n::Integer) where {F} =
    Ones{F}()

Base.similar(w::Const{F}, n::Integer) where {F} =
    Const(w.w)

Base.similar(::Vec{F}, n::Integer) where {F} =
    Vec(Vector{F}(undef, n))

Base.similar(::SVec{F}, n::Integer) where {F} =
    Vec(Array{F, 2}(undef, n, 2))

@inline is_const(::Type{SVec}) = false
@inline is_const(::Type{SVec{F}}) where {F} = false

@inline is_const(::Type{Vec}) = false
@inline is_const(::Type{Vec{F}}) where {F} = false

@inline is_const(::Type{Ones}) = true
@inline is_const(::Type{Ones{F}}) where {F} = true

@inline is_const(::Type{Const}) = true
@inline is_const(::Type{Const{F}}) where {F} = true

@inline is_const(t::Weights{F}) where {F} = is_const(typeof(t))

@inline is_const(::Type{T}) where {T<:Weights} = error("not implemented for $T)")
@inline is_const(::Type{T}) where {F, T<:Weights{F}} =
    error("not implemented for $T)")

@deprecate create_weights(x) Weights(x)

Base.size(::Ones{F}) where {F} = (0,)
Base.size(::Const{F}) where {F} = (0,)
Base.size(a::Vec{F}) where {F} = size(a.a)
Base.size(a::SVec{F}) where {F} = size(a.a)

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
"""
	Structure representing context-sensitive grammar
	Extends ExprRules.Grammar with constraints
"""
struct ContextSensitiveGrammar <: Grammar
	rules::Vector{Any}
	types::Vector{Symbol}
	isterminal::BitVector
	iseval::BitVector
	bytype::Dict{Symbol, Vector{Int}}
	childtypes::Vector{Vector{Symbol}}
	constraints::Vector{Constraint}
end

macro csgrammar(ex)
	rules = Any[]
	types = Symbol[]
	bytype = Dict{Symbol,Vector{Int}}()
	for e in ex.args
	    if isa(e, Expr)
		if e.head == :(=)
		    s = e.args[1] # name of return type
		    rule = e.args[2] # expression?
		    rvec = Any[]
		    _parse_rule!(rvec, rule)
		    for r in rvec
			push!(rules, r)
			push!(types, s)
			bytype[s] = push!(get(bytype, s, Int[]), length(rules))
		    end
		end
	    end
	end
	alltypes = collect(keys(bytype))
	is_terminal = [isterminal(rule, alltypes) for rule in rules]
	is_eval = [iseval(rule) for rule in rules]
	childtypes = [get_childtypes(rule, alltypes) for rule in rules]
	return ContextSensitiveGrammar(rules, types, is_terminal, is_eval, bytype, childtypes, [])
end

"""
    Add constraint to the grammar
"""
addconstraint!(grammar::ContextSensitiveGrammar, cons::Constraint) = push!(grammar.constraints, cons)


text:
#= Dispatch on tuples of types
dcjones commented on Apr 22, 2015

https://github.com/JuliaLang/julia/issues/10947 =# #status closed (with errors) [& multiple open issues]

#=buggyline type A end; # since 0.3! =#

 isa(A, Type{A}); #true

 isa((A,), Tuple{Type{A}}); #false

 isa(convert(Tuple{Type{A}}, (A,)), Tuple{Type{A}}); #false

#=Buggylinetype A end  #UncommentMe  =#

isa(A, Type{A}) #= true #A is not defined =#

isa((A,), (Type{A},)) # true ERROR A is not defined =#

#=To overcome the breaking change of (Normal,Normal) needing to be Tuple{Normal,Normal} I'll try to add another method that dispatches one to the other. =#
Tuple{Normal,Normal} #= ERROR: UndefVarError: Normal not defined

now known as:=#
Tuple{Number,Number} #ok
#=
ERROR: UndefVarError: Normal not defined
Stacktrace:
 [1] top-level scope
   @ REPL[30]:1=#

# Tuple{Normal,Normal} #Error # was working on  commented by @sebastiang , on Apr 22, 2015 # UndefVarError: Normal not defined
Tuple{Number,Number} #ok Tuple{Number, Number} # working!  (Q. but how to ensure they've meant the same thing ?)

#=(JeffBezanson commented on Apr 22, 2015)  Yes, it should work to add a method to convert (Normal,Normal) to Tuple{Normal,Normal}, although the resulting dispatch will probably be slow.
here's nothing wrong with using tuples. In fact now using them is a great idea, since they perform well and don't require declarations
One huge benefit is that if everybody who needs "two integers" uses Tuple{Int,Int}, then [they can share lots of code]. Declaring new types every time can be highly redundant.

[JeffBezanson closed this on Apr 24, 2015]

    toivoh commented on Apr 24, 2015
So if isa shouldn't be used to check the applicability of a method signature to an argument tuple, what should? Isn't this an inconsistency in the type system? (Making tuples not completely covariant?) I can understand if it's tricky to do something about it, but it would be good to get some more perspective about how it could be allowed to behave in the long run(?)

    JeffBezanson commented on Apr 24, 2015
                [**`Tuples are still covariant`**] =#

Tuple{Type{Int}} <: Tuple{DataType} #True  # but of course not the other way around. It's just that you can't construct a value of the left-hand type(!)
#One thing worth considering is making  typeof(T) always return Type{T} for types,
typeof(T)  # Type{T} #returns: typeof(ones)
DataType{T} #  LoadError: TypeError: in Type{...} expression, expected UnionAll, got Type{DataType}
#=ERROR: TypeError: in Type{...} expression, expected UnionAll, got Type{DataType}=#

#@toivoh commented on Apr 24, 2015: and doing reflection some other way. We could maybe even have DataType{T}. There's an awful lot of
isa(T, Type) ? Type{T} : typeof(T) #ok # typeof(ones)
#But if tuples are covariant, then:=#

isa(A, Type{A}) # true

isa((A,), Tuple{Type{A}})  # ERROR: TypeError: in Type, in parameter, expected  Type, got a value of type Matrix{Float64} ##  LoadError: TypeError: in Type, in parameter, expected Type, got a value of type Matrix{Float64} #OPEN! (since 0.4!)
#(right?) which is clearly not the case currently in 0.4, according to @dcjones original example above. =#

#= sebastiang commented on Apr 24, 2015 #  LoadError: TypeError: in Type, in parameter, expected Type, got a value of type Matrix{Float64}#I suspect this is perhaps only tangentially related, but it seems a reasonable place to ask. I'm trying to update https://github.com/timholy/HDF5.jl to use the new syntax and am running into problems I can summarize here. Given
=#
struct Y{T<:Tuple{Vararg{AbstractVector}}} end #ok now
#  Why can I not construct the type as I expect here
Y{(Vector{Int},)}() # ERROR: TypeError: in Type, in parameter, expected Type, got a value of type Tuple{DataType}

#ERROR: TypeError: apply_type: in Y, expected Type{T}, got Tuple{DataType}
Y{(Vector{Int}, Vector{String})}()

#ERROR: TypeError: apply_type: in Y, expected Type{T}, got Tuple{DataType,DataType}
#I'm hoping this is just a syntax confusion on my part.

#Same Return:
Y{Tuple{Vector{Int}}}() # Y{Tuple{Vector{Int64}}}()
Y{Tuple{Array{Int64,1}}}() # Y{Tuple{Vector{Int64}}}() # same RETURN
# But the distinction in this case between the type of the tuple and the tuple itself is confusing for me. I've done one or the other wrong.

#----------
#= sebastiang commented on Apr 24, 2015
# Tim has code which returns a tuple of dimensions. e.g. (10,5,2) for describing fixed-size arrays. Ideally, we'd want to be able to parameterize a type with a tuple of ints, e.g. FixedArray{(10,5,2)}. I can do this in v0.4 by just declaring immutable FixedArray{T}, but I don't know how to restrict the type. This is an error: =#
struct FixedArray{T<:Tuple{Vararg{Int}}} end #ok
FixedArray{(1, 2, 3)} # ERROR: TypeError: FixedArray: in T, expected T<:Tuple{Vararg{Int64}}, got Tuple{Int64,Int64,Int64}

Vararg{Int64} #} # , got Tuple{Int64,Int64,Int64} #Vararg{Int64, N} where N #now
#--
#carnaval commented on Apr 24, 2015 #What you would want is something like immutable A{T::Tuple{Vararg{Int}}} not <:. It does not work (syntax error) but probably could without too much effort. Jeff ?
#--- (JeffBezanson commented on Apr 24, 2015)

#@toivoh No, that's not what covariance means.
#@sebastiang ( )
#is now never syntax for a type. If you're writing a type, use

Tuple{}:Y{Tuple{Vector{Int},}} #ERROR: LoadError: MethodError: no method matching  -(::Type{Y{Tuple{Vector{Int64}}}}, ::Type{Tuple{}}) #OPEN!
#---
#JeffBezanson commented on Apr 24, 2015 #  `isa` -restrictions on type parameters are not implemented yet, but have come up from time to time. I was pretty sure there was an issue about it, but I can't find it. (me:lol!)

#= sebastiang commented on Apr 24, 2015 But how do I parameterize a type with a vararg tuple then?
=#
struct A{T<:Tuple{Vararg{Int}}} end #ok
#----
# This is conceptually what I want, but I see why it doesn't work: (1,2,3) is a tuple of ints, not a type which describes a tuple of ints.
A{(1, 2, 3)}()  #ERROR: LoadError: TypeError: in A, in T, expected T<:Tuple{Vararg{Int64, N} where N}, got a value of type Tuple{Int64, Int64, Int64} #OPEN!
#---

#And this works, but doesn't help me in the least. Where did my (1,2,3) go?
A{Tuple{Int}}()  #ok #(was # ERROR: TypeError: in Type{...} expression, expected UnionAll, got a value of type Matrix{Float64} )
A{Tuple{Int64}}() #ok #(was # LoadError: TypeError: in Type{...} expression, expected UnionAll, got a value of type Matrix{Float64})

#So I'm still confused as to how I can retain the parameterization by a value
(1, 2, 3)  #and restrict that value to be the right kind of Tuple, i.e. such that
isa((1, 2, 3), Tuple{Vararg{Int}}) #ok
#----Killer ERROR DETECTED
#sebastiang commented on Apr 24, 2015 @JeffBezanson, [I guess this is what your isa-restriction comment meant]; perhaps it's not possible to describe this restriction in v0.4 and it worked 'by accident' in v0.3 because tuples of types were also types? My head spins; thanks for your patience.
#Really all we want is FixedArray{T::Int...}
FixedArray{T::Int...}; #=ERROR: LoadError: UndefVarError: FixedArray not defined #solved (below):
# JeffBezanson commented on Apr 24, 2015) #That's right; we've never had isa restrictions. This doesn't work in 0.3 either:
#solution: =# #omment (for not ) #UncommentMe

#struct FixedArray{T<:(Int...)} end; #OPEN! # ERROR: LoadError: MethodError: no method matching iterate(::Type{Int64})#Reason:cannot declare constant (#reason redefinition, if not try removing constant ) #potential_ERROR

#iterate( ::T) where T <:Union{(Int..)) end; #Union
iterate(::Type{Int64}) # ::Type
#= Closest candidates are: [Julia provides solutions ]
 ```
 iterate(::Union{LinRange, StepRangeLen}) #at range.jl:664
  iterate(::Union{LinRange, StepRangeLen}, ::Int64) #at range.jl:664
  iterate(::T) where T<:Union{Base.KeySet{var"#s77", var"#s76"} ##where {var"#s77", var"#s76"<:Dict}, Base.ValueIterator{var"#s75"} where var"#s75"<:Dict} at dict.jl:693
iterate(::Type{Int64})
```
#  ERROR: syntax: invalid "::" syntax around REPL[38]:1
FixedArray{(1, 2, 3)}
#= ERROR: LoadError: TypeError: in FixedArray, in T,  expected T<:Tuple{Vararg{Int64, N} where N}, got a value of type Tuple{Int64, Int64, Int64}= #
@simonbyrne simonbyrne mentioned this issue on Dec 1, 2020
Can't dispatch on Type{} simonbyrne/KeywordDispatch.jl#10
# Open !
#----
(sebastiang commented on Apr 24, 2015)
@JeffBezanson, I guess this is what your isa-restriction comment meant; perhaps it's not possible to describe this restriction in v0.4 and it worked 'by accident' in v0.3 because tuples of types were also types? My head spins; thanks for your patience.
Really all we want is FixedArray{T::Int...}

(JeffBezanson commented on Apr 24, 2015)
That's right; we've never had isa restrictions. This doesn't work in 0.3 either:
=#
#
FixedArray{T<:(Int...)} #} #end # ERROR: UndefVarError: T not defined #Open(???)

#= sebastiang commented on Apr 24, 2015
# This is why in HDF5, the wily Tim created a sentinel type DimSize{N} to hold those integers:

# Stub types to encode fixed-size arrays for H5T_ARRAY
=#

#Buggy(keep commented): #struct  DimSize{N}; end  # Int-wrapper (can't use tuple of Int as param)
#

#struct  FixedArray{T,D<:(DimSize...)}; end =#
dimsize{N}(::Type{DimSize{N}}) = N ;
# = ERROR: LoadError: syntax: missing comma or } in argument list #readble error! =#

size{T,D}(::Type{FixedArray{T,D}}) = map(dimsize, D)::(Int...) #= ERROR: LoadError: syntax: missing comma or } in argument list =#

eltype{T,D}(::Type{FixedArray{T,D}}) = T; #= undefined T (# requires more checking )

#= toivoh commented on Apr 26, 2015
@JeffBezanson: You're right, that's not what covariance means. Thinking about it some more, I thought that was what tuple (actually 1-tuple) types mean.
Anyway, making  =# =#
typeof(T); ## = always return Type{T} for types definitely has some appeal. It would make things more consistent, e.g. since it would make

isa(x, T) isa typeof(x) <: T;

#=typeof(x) <: T =#
#=equivalent. I have no idea about the implications when it comes to performance, breakage, etc though=#

#=@JeffBezanson JeffBezanson mentioned this issue on Jun 2, 2015
Dispatch on tuples with a type #11535  #Closed !

# @cstjean cstjean mentioned this issue on Apr 4, 2017
 Segfault on tuple-type return type declaration #21271 #Closed
=#
 #still Open : (& related):
#=redesign=#

 typeof(::Type); #OPEN! 29368  StephenVavasis mentioned this issue on Oct 5, 2018
#https://github.com/JuliaLang/julia/issues/29368

#=
Can't dispatch on Type{}
https://github.com/simonbyrne/KeywordDispatch.jl/issues/10
=#
#Corrected & closed by  Simon Byrne simonbyrne #Software developer @CliMA  Pasadena, CA

(Float64,) isa Tuple{Type{Float64}}  #false

max_stars_repo_name: adamwillisMastery/MiniGtk

max_issues_repo_name: adamwillisMastery/MiniGtk

max_forks_repo_name: adamwillisMastery/MiniGtk


text:
function test_lifted_representation()
    change_fulldim_test(LiftedHRepresentation)
    change_fulldim_test(LiftedVRepresentation)

    isempty_vtest(LiftedVRepresentation(zeros(0, 2)), 0, 0)
    isempty_vtest(LiftedVRepresentation([1 0; 1 1]), 0, 2)
    isempty_vtest(LiftedVRepresentation([0 1; 0 2]), 2, 0)
    isempty_vtest(LiftedVRepresentation([0 0; 1 1]), 1, 1)

    isempty_htest(LiftedHRepresentation(ones(0, 2)), 0, 0)
    isempty_htest(LiftedHRepresentation([0 1; 0 2]), 0, 2)
    isempty_htest(LiftedHRepresentation([0 1; 0 2], BitSet(1:2)), 2, 0)
    isempty_htest(LiftedHRepresentation([0 1; 0 2], BitSet([2])), 1, 1)

    @testset "with bad arguments" begin
        A = [1 -1 -1; 0 1 0; 0 0 1]
        ls = BitSet([1])

        @test_throws ErrorException LiftedHRepresentation(A, BitSet([4]))
        ine = copy(LiftedHRepresentation(A, ls))
        @test ine.A == A
        @test ine.A !== A
        #@test linset(ine) == ls
        @test ine.linset !== ls
        @test Polyhedra.similar_type(LiftedHRepresentation{Int, Matrix{Int}}, Float64) == LiftedHRepresentation{Float64, Matrix{Float64}}
        @test Polyhedra.similar_type(LiftedHRepresentation{Int, SparseMatrixCSC{Int, Int}}, 3, Float64) == LiftedHRepresentation{Float64, SparseMatrixCSC{Float64, Int}}

        A2 = [1 1; -1 0; 0 -1]
        b2 = [1, 0, 0]
        linset2 = BitSet([1])
        ine2 = hrep(A2, b2, linset2)

        ine = LiftedHRepresentation(ine2)
        @test ine.A == A
        @test ine.linset == ls

        V = [1 0 1; 1 1 0]
        Vlinset = BitSet(2)
        @test_throws ErrorException LiftedVRepresentation(V, BitSet([4]))
        ext = copy(LiftedVRepresentation(V, Vlinset))
        @test ext.R == V
        @test ext.R !== V
        #@test linset(ext) == Vlinset
        @test ext.linset !== Vlinset

        @test Polyhedra.similar_type(LiftedVRepresentation{Int, SparseMatrixCSC{Int, Int}}, Float64) == LiftedVRepresentation{Float64, SparseMatrixCSC{Float64, Int}}
        @test Polyhedra.similar_type(LiftedVRepresentation{Int, Matrix{Int}}, 3, Float64) == LiftedVRepresentation{Float64, Matrix{Float64}}
    end
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
"""
Propagation module
"""

# struct for handling information of the problem
struct ShootingSettings
	n::Int
	nr::Int
	n_sv::Int
	prob_stm::ODEProblem
	tofs::Array
	r0::Array    # fixed vector
	rf::Array    # fixed vector
	v0::Array
	vf::Array
	method
	reltol::Float64
	abstol::Float64
	tolDC::Float64
	tolConv::Float64
	maxiter::Int
	verbosity::Int
end

"""
Two-stage shooting algorithm for n-impulsive trajectory design
"""
function twostage_shooting(n_sv::Int, svs::Array, tofs::Array, prob_stm::ODEProblem; kwargs...)

	# FIXME -- take from kwargs
	reltol = 1.e-12
	abstol = 1.e-12
	method = Tsit5()
	tolDC = 1.e-11
	tolConv = 1.e-6
	maxiter = 15
	verbosity = 1

	n  = length(svs)  # number of nodes
	nr = n_sv ÷ 2     # number of elements of positions (== length of velocities)

	println("Using $n nodes")
	for (i,sv) in enumerate(svs)
		print("Node $i: \n")
		println(sv)
	end

	# initial and final position and velocity vectors
	r0 = svs[1][1:nr]
	v0 = svs[1][nr+1:end]
	rf = svs[end][1:nr]
	vf = svs[end][nr+1:end]

	# inner-loop decision vector: velocities at first n-1 nodes
	x_inner = zeros(nr*(n-1))
	for i = 1:n-1
		x_inner[1+(i-1)*nr:i*nr] = svs[i][nr+1:end]
	end

	# outer-loop decision vector
	x_outer = zeros((n-2)*nr)
	for j = 1:n-2
		x_outer[1+(j-1)*nr:j*nr] = svs[j+1][1:nr]
	end

	# construct ShootingSettings
	Settings = ShootingSettings(n, nr, n_sv, prob_stm, tofs, r0, rf, v0, vf,
		method, reltol, abstol, tolDC, tolConv, maxiter, verbosity)

	# initialize storage
	J_inner = zeros(nr*(n-1), nr*(n-1))

	# outer-loop shooting to minimize velocity discontinuity
	outer_loop_shooting!(x_outer, x_inner, J_inner, Settings)

	# re-construct nodes of final solution
	conv_rs = vcat(r0, x_outer, rf)
	conv_vs = vcat(x_inner, Settings.vf)

	nodes_conv = []
	for i = 1:n
		sv = vcat( conv_rs[1+(i-1)*nr : i*nr] , conv_vs[1+(i-1)*nr : i*nr] )
		push!(nodes_conv, sv)
	end
	return nodes_conv
end

"""
Outer-loop finds least-square solution to minimize velocity discontinuity.
Function mutates `x0_outer`
"""
function outer_loop_shooting!(x0_outer, x0_inner, J_inner, Settings::ShootingSettings, verbose::Bool=true)

	# initialize old velocities
	old_vs = vcat(x0_inner, Settings.vf)
	old_cost = 0.0

	for i_outer = 1:Settings.maxiter

		# inner-loop shooting to ensure position continuity
		inner_loop_shooting!(x0_outer, x0_inner, J_inner, Settings, true)

		# compute new velocities
		new_vs = vcat(x0_inner, Settings.vf)

		# break if norm of cost stops improving
		dvs = new_vs - old_vs
		cost = norm(dvs)

		if verbose==true
			println("Outer iteration $i_outer ... cost: $cost")
		end

		if (i_outer > 1) && (old_cost - cost < Settings.tolConv)
			println("Outer-loop achieved tolerance!")
			break
		end

		# compute J_outer via finite difference
		function g(x_outer)
			inner_loop_shooting!(x_outer, x0_inner, J_inner, Settings, false)
			new_vs = vcat(x0_inner, Settings.vf)
			return new_vs - old_vs
		end
		J_outer = FiniteDiff.finite_difference_jacobian(g, x0_outer)

		# least-square update
		# print("x0_outer: "); println(length(x0_outer))
		# print("J_outer: "); println(size(J_outer))
		# print("new_vs: "); println(length(new_vs))

		# updates
		x0_outer[:] = x0_outer -inv(transpose(J_outer)*J_outer) * transpose(J_outer) * new_vs
		old_vs[:] = new_vs
		old_cost = cost
	end

	return
end

"""
Inner-loop corrects velocity vectors to ensure position continuity.
Function mutates `x0_inner`
"""
function inner_loop_shooting!(x0_outer, x0_inner, J_inner, Settings::ShootingSettings, verbose::Bool=false)

	r_i2_prop = zeros(Settings.nr * (Settings.n-1))

	for i_inner = 1:Settings.maxiter

		# propagate n nodes
		for i = 1:Settings.n-1
			# re-construct array of positions from node 1 ~ n-1
			rs_vec = vcat(Settings.r0, x0_outer)
			#
			r_i = rs_vec[1+(i-1)*Settings.nr : i*Settings.nr]
			v_i = x0_inner[1+(i-1)*Settings.nr : i*Settings.nr]
			x0i =  vcat(r_i, v_i, reshape(I(Settings.n_sv), (Settings.n_sv*Settings.n_sv)))[:]
			_prob = remake(Settings.prob_stm; tspan=(0.0, Settings.tofs[i]), u0=x0i)
	        sol = DifferentialEquations.solve(_prob, Settings.method, reltol=Settings.reltol, abstol=Settings.abstol)

	        # store results
	        r_i2_prop[1+(i-1)*Settings.nr : i*Settings.nr] = sol.u[end][1:Settings.nr]

	    	# fill-in upper-left submatrix of STM into Jacobian
	    	STM = transpose( reshape(sol.u[end][Settings.n_sv+1:end], (Settings.n_sv, Settings.n_sv)) )
	    	J_inner[1+(i-1)*Settings.nr:i*Settings.nr, 1+(i-1)*Settings.nr:i*Settings.nr] = -STM[1:Settings.nr, Settings.nr+1:end]  # FIXME
		end

		# compute final position offset: δr_i2 = r_i2_guess - r_i2_prop
		δr_i2 = vcat(x0_outer, Settings.rf) - r_i2_prop

		# check breaking condition
		err = norm(δr_i2)
		if verbose==true
			println("   Inner iteration $i_inner ... err: $err")
		end
		if err < Settings.tolDC
			if verbose==true
				println("   Inner-loop achieved tolerance!")
			end
			break
		end

		# correct velocity: v_i = v_i_guess - J_inner^-1 * δr_i2
		x0_inner[:] = x0_inner - inv(J_inner) * δr_i2

	end

	return
end


text:
"""
    function randomAdapt_(p, code, platform)

Function that adapt the random function invocations of the code to be executable in the different platforms.
"""
function randomAdapt_(p::Program_, code::Expr, platform::String)

    if platform == "cpu"
        for i in VALIDDISTRIBUTIONS
            code = postwalk(x -> @capture(x,$i(v__)) ? :(rand($i($(v...)))) : x, code)
        end
    elseif platform == "gpu"
        for i in VALIDDISTRIBUTIONS

            if inexpr(code,i) && !(i in VALIDDISTRIBUTIONSCUDA)
                error(i," random distribution valid for cpu but still not implemented in gpu.")
            else

                s = Meta.parse(string("AgentBasedModels.",i,"CUDA"))
                code = postwalk(x -> @capture(x,$i(v__)) ? :($s(rand(),$(v...))) : x, code)

            end

        end
    end

    return code
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
import Calculus:differentiate,jacobian,simplify,SymbolParameter,isminus

function simplify(::SymbolParameter{:-}, args)
	if length(args) == 2 && args[2]==0
		return args[1]
	end
	# Remove any 0's in a subtraction
    args = map(simplify, filter(x -> x != 0, args))
    if length(args) == 0
        return 0
    # Special Case: simplify(:(x - x)) == 0
    elseif length(args) == 2 && args[1] == args[2]
        return 0
    # Special Case: simplify(:(x - (-y))) == x + y
    elseif length(args) == 2 && isminus(args[2])
        return Expr(:call, :+, args[1], args[2].args[2])
    else
	        return Expr(:call, :-, args...)
    end
end

function differentiate(ex::Expr,wrt::Expr)
	if ex.head==:vect || ex.head ==:vcat
		return differentiate(SymbolParameter(:vect), ex.args[1:end], wrt)
	elseif ex.head == :ref
        return ex==wrt ? 1 : 0
    elseif ex.head != :call
		error("Unrecognized expression $ex")
	end
    return simplify(differentiate(SymbolParameter(ex.args[1]), ex.args[2:end], wrt))
end

function differentiate(ex::Expr, targets::Vector{Expr})
	n = length(targets)
	exprs = Array(Any, n)
	for i in 1:n
		exprs[i] = differentiate(ex, targets[i])
	end
	return exprs
end

differentiate(ex::Union{Number,Symbol},wrt::Expr) = 0

function differentiate(::SymbolParameter{:vect}, args, wrt)
	for i = 1:length(args)
		args[i] = differentiate(args[i],wrt)
	end
	return Expr(:vect,args...)
end

function jacobian(ex::Expr, targets::Vector)
    @assert ex.head==:vect || ex.head==:vcat
    exprs = Expr(:vcat)
    for i = 1:length(ex.args)
        push!(exprs.args,Expr(:row))
        for j = 1:length(targets)
            push!(exprs.args[i].args,simplify(differentiate(ex.args[i],targets[j])))
        end
    end
    return exprs
end

differentiate(::SymbolParameter{:max}, args, wrt) = Expr(:if,:($(args[1])>$(args[2])),differentiate(args[1],wrt),differentiate(args[2],wrt))
differentiate(::SymbolParameter{:min}, args, wrt) = Expr(:if,:($(args[1])<$(args[2])),differentiate(args[1],wrt),differentiate(args[2],wrt))

function simplify(ex::Expr)
	if in(ex.head,[:vcat,:row])
		for i = 1:length(ex.args)
			ex.args[i]=simplify(ex.args[i])
		end
		return ex
	end

	if ex.head != :call
        return ex
    end

    if all(Calculus.isnumber, ex.args[2:end]) && length(ex.args) > 1
        return eval(current_module(), ex)
    end
    new_ex = simplify(SymbolParameter(ex.args[1]), ex.args[2:end])
    while new_ex != ex
        new_ex, ex = simplify(new_ex), new_ex
    end
    return new_ex
end


text:
A = randn(3,5)
b = randn(3)

x_pg = Variable(5)
prob_pg = problem(ls(A*x_pg - b) + 1e-3*norm(x_pg, 1))
sol_pg = solve(prob_pg, StructuredOptimization.PG())

x_zfpr = Variable(5)
prob_zfpr = problem(ls(A*x_zfpr - b) + 1e-3*norm(x_zfpr, 1))
sol_zfpr = solve(prob_zfpr, StructuredOptimization.ZeroFPR())

x_pnc = Variable(5)
prob_pnc = problem(ls(A*x_pnc - b) + 1e-3*norm(x_pnc, 1))
sol_pnc = solve(prob_pnc, StructuredOptimization.PANOC())


text:
using ProxBTF
σ   = 3# noise s.d.
x   = [1:1.:100;]
fx  = 13 .* sin.((4*π/100).*x)
y = fx + σ .* randn(100)
result_summary,result_quantile,problem,elapse_time=pbtf(y,x,2,nsample=3000)
visualize(result_quantile,problem,legend_position=:top)
σ   = 1.0# noise s.d.
x   = [0.1:0.1:10;]
fx  = x .+ sin.(x)
y = fx + σ .* randn(100)
result_summary,result_quantile,problem,elapse_time=pbsrtf(y,x,2,"increasing")
visualize(result_quantile,problem)
fx = map(xi -> 0<= xi <=2 ? 10-5xi : 2<xi <= 8 ? 0 : 5xi-40,x)
y = fx + σ .* randn(100)
result_summary,result_quantile,problem,elapse_time=pbsrtf(y,x,1,"convex")
visualize(result_quantile,problem)
fx = map(xi -> 0<= xi <= 8 ? 0 : (xi-8)^3,x)
y = fx + σ .* randn(100)
result_summary,result_quantile,problem,elapse_time=pbsrtf(y,x,2,"inc-convex")
visualize(result_quantile,problem)


text:
using Base.Test
using Base.Iterators

@testset "Weber" begin
  @testset "Trial Sequencing" begin
    include("test_moment_timing.jl")
    include("test_moment_sequencing.jl")
    include("test_moment_indexing.jl")
    include("test_compound_moments.jl")
    include("test_moment_preparation.jl")
    include("test_moment_looping.jl")
    include("test_moment_conditions.jl")
  end
  include("test_record_columns.jl")
  include("test_moment_checks.jl")
  include("test_extensions.jl")
  include("test_oddball.jl")
end


text:
import Base: getindex, setindex!, axes

axes(x::TS) = axes(x.values)
axes(x::TS, d) = axes(x.values, d)

# all element indexing
getindex(x::TS) = x
getindex(x::TS, ::Colon) = x
getindex(x::TS, ::Colon, ::Colon) = x

# cartesian index support
getindex(x::TS, idx::CartesianIndex{2}) = x.values[idx.I[1], idx.I[2]]
getindex(x::TS, idx::CI) where {CI<:AbstractArray{CartesianIndex{2}}} = TS(x.values[idx],
                                                                           x.index[[idx[i].I[1] for i in 1:length(idx)]],
                                                                           x.fields[unique([idx[i].I[2] for i in 1:length(idx)])])

# general interface
getindex(x::TS, r, c) = (rows=findrows(r,x); cols=findcols(c,x); TS(x.values[rows,cols],x.index[rows],x.fields[cols]))
getindex(x::TS, r, ::Colon) = (rows=findrows(r,x); TS(x.values[rows,:],x.index[rows],x.fields))
getindex(x::TS, ::Colon, c) = (cols=findcols(c,x); TS(x.values[:,cols],x.index,x.fields[cols]))

# types only used to index rows
getindex(x::TS, r::R) where {R<:Int} = (rows=findrows(r,x); TS(x.values[rows,:],x.index[rows],x.fields))
getindex(x::TS, r::R) where {R<:AbstractVector{<:Integer}} = (rows=findrows(r,x); TS(x.values[rows,:],x.index[rows],x.fields))
getindex(x::TS, r::R) where {R<:TimeType} = (rows=findrows(r,x); TS(x.values[rows,:],x.index[rows],x.fields))
getindex(x::TS, r::R) where {R<:AbstractVector{<:TimeType}} = (rows=findrows(r,x); TS(x.values[rows,:],x.index[rows],x.fields))
getindex(x::TS, r::R) where {R<:AbstractString} = (rows=findrows(r,x); TS(x.values[rows,:],x.index[rows],x.fields))

# types only used to index columns
getindex(x::TS, c::C) where {C<:Symbol} = (cols=findcols(c,x); TS(x.values[:,cols],x.index,x.fields[cols]))
getindex(x::TS, c::C) where {C<:AbstractVector{<:Symbol}} = (cols=findcols(c,x); TS(x.values[:,cols],x.index,x.fields[cols]))

setindex!(x::TS, v, ::Colon) = x.values[:] .= v
setindex!(x::TS, v, ::Colon, ::Colon) = x[:] = v
setindex!(x::TS, v, ::Colon, c) = x.values[:, findcols(c,x)] .= v
setindex!(x::TS, v, r, ::Colon) = x.values[findrows(r,x), :] .= v
setindex!(x::TS, v, r, c) = x.values[findrows(r,x), findcols(c,x)] .= v

setindex!(x::TS, v, r::R) where {R<:TimeType} = x.values[findrows(r,x), :] .= v
setindex!(x::TS, v, r::R) where {R<:AbstractVector{<:TimeType}} = x.values[findrows(r,x), :] .= v
setindex!(x::TS, v, r::R) where {R<:AbstractString} = x.values[findrows(r,x), :] .= v
setindex!(x::TS, v, r::R) where {R<:Int} = x.values[findrows(r,x), :] .= v
setindex!(x::TS, v, r::R) where {R<:AbstractVector{<:Integer}} = x.values[findrows(r,x), :] .= v

# single column
makecolumns(x::TS, v::V, c::C) where {V<:VALARR, C<:FLDTYPE} = TS(v, x.index, [c])
makecolumns(x::TS, v::V, c::C) where {V<:TS, C<:FLDTYPE} = TS(v.values, x.index, [c])
makecolumns(x::TS, v::V, c::C) where {V<:VALTYPE, C<:FLDTYPE} = TS(zeros(typeof(v), size(x,1)).+v, x.index, [c])
# multiple columns
makecolumns(x::TS, v::V, c::C) where {V<:VALARR, C<:FLDARR} = TS(v, x.index, c)
makecolumns(x::TS, v::V, c::C) where {V<:TS, C<:FLDARR} = TS(v.values, x.index, c)

# column modification/addition via e.g. X[:Column] = randn(N)
function setindex!(x::TS, v::V, c::C) where {V<:Union{<:VALARR,VALTYPE,TS},C<:Union{AbstractVector{Symbol},Symbol}}
    found = sum(findcols(c,x)) > 0
    if found
        x.values[:, findcols(c,x)] .= v
    else
        y = [x makecolumns(x, v, c)]
        x.values = y.values
        x.fields = y.fields
    end
end


text:
using Kinetic
using KitBase.ProgressMeter

cd(@__DIR__)
ks, ctr, face, simTime = Kinetic.initialize("briowu_1d.txt")

dt = Kinetic.timestep(ks, ctr, simTime)
nt = Int(floor(ks.set.maxTime / dt))+1
res = zeros(5, 2)

KS = ks
@showprogress for iter in 1:nt
    #dt = timestep(KS, ctr, simTime)
    Kinetic.reconstruct!(ks, ctr)
    Kinetic.evolve!(ks, ctr, face, dt; mode=:kcu, isPlasma=true)
    #=
    @inbounds Threads.@threads for i = 1:KS.pSpace.nx+1
        flux_kfvs!(
            face[i].fw,
            face[i].fh0,
            face[i].fh1,
            face[i].fh2,
            face[i].fh3,
            ctr[i-1].h0 .+ 0.5 .* ctr[i-1].dx .* ctr[i-1].sh0,
            ctr[i-1].h1 .+ 0.5 .* ctr[i-1].dx .* ctr[i-1].sh1,
            ctr[i-1].h2 .+ 0.5 .* ctr[i-1].dx .* ctr[i-1].sh2,
            ctr[i-1].h3 .+ 0.5 .* ctr[i-1].dx .* ctr[i-1].sh3,
            ctr[i].h0 .- 0.5 .* ctr[i].dx .* ctr[i].sh0,
            ctr[i].h1 .- 0.5 .* ctr[i].dx .* ctr[i].sh1,
            ctr[i].h2 .- 0.5 .* ctr[i].dx .* ctr[i].sh2,
            ctr[i].h3 .- 0.5 .* ctr[i].dx .* ctr[i].sh3,
            KS.vSpace.u,
            KS.vSpace.weights,
            dt,
            ctr[i-1].sh0,
            ctr[i-1].sh1,
            ctr[i-1].sh2,
            ctr[i-1].sh3,
            ctr[i].sh0,
            ctr[i].sh1,
            ctr[i].sh2,
            ctr[i].sh3,
        )
    end
    @inbounds Threads.@threads for i = 1:KS.pSpace.nx+1
        flux_em!(
            face[i].femL,
            face[i].femR,
            ctr[i-2].E,
            ctr[i-2].B,
            ctr[i-1].E,
            ctr[i-1].B,
            ctr[i].E,
            ctr[i].B,
            ctr[i+1].E,
            ctr[i+1].B,
            ctr[i-1].ϕ,
            ctr[i].ϕ,
            ctr[i-1].ψ,
            ctr[i].ψ,
            ctr[i-1].dx,
            ctr[i].dx,
            KS.gas.Ap,
            KS.gas.An,
            KS.gas.D,
            KS.gas.sol,
            KS.gas.χ,
            KS.gas.ν,
            dt,
        )
    end
    =#
    #Kinetic.update!(ks, ctr, face, dt, res; coll=:bgk, bc=:extra, isMHD=true)
    # it's equivalent to the following process
    sumRes = zeros(5, 2)
    sumAvg = zeros(5, 2)
    Threads.@threads for i = 1:ks.pSpace.nx
        Kinetic.step!(
            ks,
            face[i],
            ctr[i],
            face[i+1],
            dt,
            sumRes,
            sumAvg,
        )
    end

    for i = 1:2
        ctr[1-i].w .= ctr[1].w
        ctr[1-i].prim .= ctr[1].prim
        ctr[1-i].h0 .= ctr[1].h0
        ctr[1-i].h1 .= ctr[1].h1
        ctr[1-i].h2 .= ctr[1].h2
        ctr[1-i].E .= ctr[1].E
        ctr[1-i].B .= ctr[1].B
        ctr[1-i].ϕ = deepcopy(ctr[1].ϕ)
        ctr[1-i].ψ = deepcopy(ctr[1].ψ)
        ctr[1-i].lorenz .= ctr[1].lorenz

        ctr[KS.pSpace.nx+i].w .= ctr[KS.pSpace.nx].w
        ctr[KS.pSpace.nx+i].prim .= ctr[KS.pSpace.nx].prim
        ctr[KS.pSpace.nx+i].h0 .= ctr[KS.pSpace.nx].h0
        ctr[KS.pSpace.nx+i].h1 .= ctr[KS.pSpace.nx].h1
        ctr[KS.pSpace.nx+i].h2 .= ctr[KS.pSpace.nx].h2
        ctr[KS.pSpace.nx+i].E .= ctr[KS.pSpace.nx].E
        ctr[KS.pSpace.nx+i].B .= ctr[KS.pSpace.nx].B
        ctr[KS.pSpace.nx+i].ϕ = deepcopy(ctr[KS.pSpace.nx].ϕ)
        ctr[KS.pSpace.nx+i].ψ = deepcopy(ctr[KS.pSpace.nx].ψ)
        ctr[KS.pSpace.nx+i].lorenz .= ctr[KS.pSpace.nx].lorenz
    end
end

sol = zeros(ks.pSpace.nx, 10, 2)
for i in 1:ks.pSpace.nx
    sol[i, 1, 1] = ctr[i].prim[1,1]
    sol[i, 1, 2] = ctr[i].prim[1,2] / ks.gas.me
    sol[i, 2:4, 1] .= ctr[i].prim[2:4,1]
    sol[i, 2:4, 2] .= ctr[i].prim[2:4,2]
    sol[i, 5, 1] = 1. / ctr[i].prim[5,1]
    sol[i, 5, 2] = ks.gas.me / ctr[i].prim[5,2]

    sol[i, 6, 1] = ctr[i].B[2]
    sol[i, 6, 2] = ctr[i].E[1]
end
using Plots
plot(ks.pSpace.x[1:ks.pSpace.nx], sol[:,1,:])
plot(ks.pSpace.x[1:ks.pSpace.nx], sol[:,6,:])

using JLD2
@save "refsol.jld2" ks ctr


text:
m1 = fit(LinearModel, form, data)
λlm = dispersion(m1)

loss1 = RobustModels.L2Loss()
loss2 = RobustModels.TukeyLoss()
est1 = MEstimator(loss1)
est2 = MEstimator(loss2)

funcs = ( dof,
          dof_residual,
          confint,
          deviance,
          nulldeviance,
          loglikelihood,
          nullloglikelihood,
          dispersion,
          nobs,
          stderror,
          vcov,
          residuals,
          predict,
          response,
          weights,
          workingweights,
          fitted,
          predict,
          isfitted,
          islinear,
          leverage,
          modelmatrix,
          projectionmatrix,
          Estimator,
          scale
        )

@testset "linear: L2 estimator" begin
    println("\n\t\u25CF Estimator: L2")

    # OLS
    println(m1)
    println(" lm              : ", coef(m1))

    # Formula, dense and sparse entry  and methods :cg and :chol
    @testset "(type, method): ($(typeof(A)),\t$(method))" for (A, b) in ((form, data), (X, y), (sX, y)), method in (:cg, :chol)
        name  = if A==form; "formula" elseif A==X; "dense  " else "sparse " end
        name *= if method==:cg; ",  cg" else ",chol" end
        # use the dispersion from GLM to ensure that the loglikelihood is correct
        m = fit(RobustLinearModel, A, b, est1; method=method, verbose=false, initial_scale=λlm)
        β = copy(coef(m))
        println("rlm($name): ", β)
        @test_nowarn println(m)
        @test isapprox(coef(m1), β; rtol=1e-5)

        # refit
        refit!(m, y; verbose=false, initial_scale=:extrema)
        println("$m")
        @test all(coef(m) .== β)

        # make sure the methods for RobustLinearModel are well defined
        @testset "method: $(f)" for f in funcs
            if f in (Estimator, weights, workingweights, islinear, isfitted,
                     leverage, modelmatrix, projectionmatrix, scale)
                # method is not defined in GLM
                @test_nowarn f(m)
            else
                var = f(m1)
                robvar = f(m)
                if isa(var, Union{AbstractArray, Tuple})
                    if f != vcov
                        @test isapprox(var, robvar; rtol=1e-4)
                    end
                else
                    if f == dof
                        @test var == robvar + 1
                    elseif f == dof_residual
                        @test var == robvar
                    elseif f in (deviance, nulldeviance)
                        s = scale(m)
                        ## the deviance in RobustModels is the scaled deviance (divided by
                        ## the squared dispersion), whereas in GLM it is not.
                        @test isapprox(var, robvar * s^2; rtol=1e-4)
                    elseif f in (loglikelihood, nullloglikelihood)
                        ## TODO: should work
                        @test_broken isapprox(var, robvar; rtol=1e-4)
#                        @test isapprox(var + log(λ), (robvar + RobustModels.fullloglikelihood(r)) * s^2/λ^2 - log(RobustModels.estimator_norm(r.est)); rtol=1e-4)
#                        @test isapprox(var + log(dispersion(m1)), robvar + log(s); rtol=1e-4)
                    else
                        @test isapprox(var, robvar; rtol=1e-4)
                    end
                end
            end
        end

        # TauEstimator interface
        m3 = fit(RobustLinearModel, A, b, TauEstimator{TukeyLoss}(); method=method, initial_scale=λlm)
        @test_nowarn tauscale(m3)

        # later fit!
        m2 = fit(RobustLinearModel, A, b, est1; method=method, dofit=false)
        @test all(0 .== coef(m2))
        fit!(m2; verbose=false, initial_scale=:mad)
        @test all(β .== coef(m2))
    end
end


text:
"""
  module Multroot

Calculate mutiple roots of a real or complex polymomial exploiting the
multiplicity structure of the root.
Using methods develped by Zhang et. al.
"""
module Multroot

using Printf
using LinearAlgebra

export multroot, gcdroot, pejroot
export sylves, sylves1, sylmat, cauchymt, scalsq
export hqrt, hessqr, forsub, backsub
export zminsv, zminsv1
export polyzeros, polymult, polytransform

include("polyzeros.jl")

include("sylves.jl")
include("sylmat.jl")
include("polymult.jl")
include("zminsv.jl")
include("hessqr.jl")
include("hqrt.jl")
include("cauchymt.jl")
include("scalsq.jl")
include("forsub.jl")
include("backsub.jl")
include("gcdgn.jl")
include("multroot.jl")
include("pejroot.jl")
include("gcdroot.jl")

include("polyscale.jl")
include("horner.jl")

end # module


text:
mutable struct NormalizationConstant{T, PP<:AbstractQuasiMatrix{T}} <: AbstractCachedVector{T}
    P::PP # OPs
    data::Vector{T}
    datasize::Tuple{Int}

    function NormalizationConstant{T, PP}(P::PP) where {T,PP<:AbstractQuasiMatrix{T}}
        μ = inv(sqrt(sum(orthogonalityweight(P))))
        new{T, PP}(P, [μ], (1,))
    end
end

NormalizationConstant(P::AbstractQuasiMatrix{T}) where T = NormalizationConstant{T,typeof(P)}(P)

size(K::NormalizationConstant) = (∞,)

# How we populate the data
# function _normalizationconstant_fill_data!(K::NormalizationConstant, J::Union{BandedMatrix,Symmetric{<:Any,BandedMatrix},Tridiagonal,SymTridiagonal}, inds)
#     dl, _, du = bands(J)
#     @inbounds for k in inds
#         K.data[k] = sqrt(du[k-1]/dl[k]) * K.data[k-1]
#     end
# end

function _normalizationconstant_fill_data!(K::NormalizationConstant, J, inds)
    @inbounds for k in inds
        K.data[k] = sqrt(J[k,k-1]/J[k-1,k]) * K.data[k-1]
    end
end

LazyArrays.cache_filldata!(K::NormalizationConstant, inds) = _normalizationconstant_fill_data!(K, jacobimatrix(K.P), inds)

struct Normalized{T, OPs<:AbstractQuasiMatrix{T}, NL} <: OrthogonalPolynomial{T}
    P::OPs
    scaling::NL # Q = P * Diagonal(scaling)
end

normalizationconstant(P) = NormalizationConstant(P)
Normalized(P::AbstractQuasiMatrix{T}) where T = Normalized(P, normalizationconstant(P))
Normalized(Q::Normalized) = Q

struct NormalizedBasisLayout{LAY<:AbstractBasisLayout} <: AbstractBasisLayout end

MemoryLayout(::Type{<:Normalized{<:Any, OPs}}) where OPs = NormalizedBasisLayout{typeof(MemoryLayout(OPs))}()

struct QuasiQR{T, QQ, RR} <: Factorization{T}
    Q::QQ
    R::RR
end

QuasiQR(Q::AbstractQuasiMatrix{T}, R::AbstractMatrix{V}) where {T,V} =
    QuasiQR{promote_type(T,V),typeof(Q),typeof(R)}(Q, R)

Base.iterate(S::QuasiQR) = (S.Q, Val(:R))
Base.iterate(S::QuasiQR, ::Val{:R}) = (S.R, Val(:done))
Base.iterate(S::QuasiQR, ::Val{:done}) = nothing

axes(Q::Normalized) = axes(Q.P)
==(A::Normalized, B::Normalized) = A.P == B.P

# There is no point in a Normalized OP thats ==, so just return false
==(A::Normalized, B::OrthogonalPolynomial) = false
==(A::OrthogonalPolynomial, B::Normalized) = false
==(A::Normalized, B::AbstractQuasiMatrix) = false
==(A::AbstractQuasiMatrix, B::Normalized) = false
==(A::Normalized, B::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}) = false
==(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}, B::Normalized) = false

_p0(Q::Normalized) = Q.scaling[1]

# p_{n+1} = (A_n * x + B_n) * p_n - C_n * p_{n-1}
# q_{n+1}/h[n+1] = (A_n * x + B_n) * q_n/h[n] - C_n * p_{n-1}/h[n-1]
# q_{n+1} = (h[n+1]/h[n] * A_n * x + h[n+1]/h[n] * B_n) * q_n - h[n+1]/h[n-1] * C_n * p_{n-1}

function recurrencecoefficients(Q::Normalized)
    A,B,C = recurrencecoefficients(Q.P)
    h = Q.scaling
    h[2:∞] ./ h .* A, h[2:∞] ./ h .* B, Vcat(zero(eltype(Q)), h[3:∞] ./ h .* C[2:∞])
end

# x * p[n] = c[n-1] * p[n-1] + a[n] * p[n] + b[n] * p[n+1]
# x * q[n]/h[n] = c[n-1] * q[n-1]/h[n-1] + a[n] * q[n]/h[n] + b[n] * q[n+1]/h[n+1]
# x * q[n+1] = c[n-1] * h[n]/h[n-1] * q[n-1] + a[n] * q[n] + b[n] * h[n]/h[n+1] * q[n+1]

# q_{n+1}/h[n+1] = (A_n * x + B_n) * q_n/h[n] - C_n * p_{n-1}/h[n-1]
# q_{n+1} = (h[n+1]/h[n] * A_n * x + h[n+1]/h[n] * B_n) * q_n - h[n+1]/h[n-1] * C_n * p_{n-1}
function jacobimatrix(Q::Normalized)
    X = jacobimatrix(Q.P)
    a,b = X[band(0)], X[band(-1)]
    h = Q.scaling
    Symmetric(_BandedMatrix(Vcat(a', (b .* h ./ h[2:end])'), ∞, 1, 0), :L)
end

orthogonalityweight(Q::Normalized) = orthogonalityweight(Q.P)
singularities(Q::Normalized) = singularities(Q.P)

function demap(Q::Normalized)
    P,D =  arguments(ApplyLayout{typeof(*)}(), Q)
    demap(P) * D
end

# Sometimes we want to expand out, sometimes we don't

QuasiArrays.ApplyQuasiArray(Q::Normalized) = ApplyQuasiArray(*, arguments(ApplyLayout{typeof(*)}(), Q)...)

ArrayLayouts.mul(Q::Normalized, C::AbstractArray) = ApplyQuasiArray(*, Q, C)

grid(Q::SubQuasiArray{<:Any,2,<:Normalized}) = grid(view(parent(Q).P, parentindices(Q)...))

# transform_ldiv(Q::Normalized, C::AbstractQuasiArray) = Q.scaling .\ (Q.P \ C)
function transform_ldiv(Q::Normalized, C::AbstractQuasiArray)
    c = paddeddata(Q.P \ C)
    [Q.scaling[axes(c,1)] .\ c; zeros(eltype(c), ∞)]
end

function transform_ldiv(V::SubQuasiArray{<:Any,2,<:Normalized}, C::AbstractQuasiArray)
    Q = parent(V)
    P = Q.P
    kr, jr = parentindices(V)
    c = transform_ldiv(view(P, kr, jr), C)
    Q.scaling[axes(c,1)] .\ c
end

arguments(::ApplyLayout{typeof(*)}, Q::Normalized) = Q.P, Diagonal(Q.scaling)
_mul_arguments(Q::Normalized) = arguments(ApplyLayout{typeof(*)}(), Q)
_mul_arguments(Q::QuasiAdjoint{<:Any,<:Normalized}) = arguments(ApplyLayout{typeof(*)}(), Q)

# table stable identity if A.P == B.P
@inline _normalized_ldiv(An, C, Bn) = An \ (C * Bn)
@inline _normalized_ldiv(An, C::Eye{T}, Bn) where T = FillArrays.SquareEye{promote_type(eltype(An),T,eltype(Bn))}(∞)
@inline copy(L::Ldiv{<:NormalizedBasisLayout,<:NormalizedBasisLayout}) = _normalized_ldiv(Diagonal(L.A.scaling), L.A.P \ L.B.P, Diagonal(L.B.scaling))
@inline copy(L::Ldiv{Lay,<:NormalizedBasisLayout}) where Lay = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A, L.B))
@inline copy(L::Ldiv{<:NormalizedBasisLayout,Lay}) where Lay = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
@inline copy(L::Ldiv{<:NormalizedBasisLayout,Lay,<:Any,<:AbstractQuasiVector}) where Lay = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
@inline copy(L::Ldiv{Lay,<:NormalizedBasisLayout}) where Lay<:AbstractBasisLayout = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A, L.B))
@inline copy(L::Ldiv{<:NormalizedBasisLayout,Lay}) where Lay<:AbstractBasisLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
@inline copy(L::Ldiv{Lay,<:NormalizedBasisLayout}) where Lay<:AbstractLazyLayout = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A, L.B))
@inline copy(L::Ldiv{<:NormalizedBasisLayout,Lay}) where Lay<:AbstractLazyLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
@inline copy(L::Ldiv{<:NormalizedBasisLayout,Lay,<:Any,<:AbstractQuasiVector}) where Lay<:AbstractLazyLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))
@inline copy(L::Ldiv{ApplyLayout{typeof(*)},<:NormalizedBasisLayout}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A, L.B))
for Lay in (:(ApplyLayout{typeof(*)}),:(BroadcastLayout{typeof(+)}),:(BroadcastLayout{typeof(-)}))
    @eval begin
        @inline copy(L::Ldiv{<:NormalizedBasisLayout,$Lay}) = copy(Ldiv{ApplyLayout{typeof(*)},$Lay}(L.A, L.B))
        @inline copy(L::Ldiv{<:NormalizedBasisLayout,$Lay,<:Any,<:AbstractQuasiVector}) = copy(Ldiv{ApplyLayout{typeof(*)},$Lay}(L.A, L.B))
    end
end

# want to use special re-expansion routines without expanding Normalized basis
@inline copy(L::Ldiv{<:NormalizedBasisLayout,BroadcastLayout{typeof(*)}}) = copy(Ldiv{BasisLayout,BroadcastLayout{typeof(*)}}(L.A, L.B))
@inline copy(L::Ldiv{<:NormalizedBasisLayout,BroadcastLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = copy(Ldiv{BasisLayout,BroadcastLayout{typeof(*)}}(L.A, L.B))

###
# show
###
Base.array_summary(io::IO, C::NormalizationConstant{T}, inds) where T = print(io, "NormalizationConstant{$T}")
show(io::IO, Q::Normalized) = print(io, "Normalized($(Q.P))")
show(io::IO, ::MIME"text/plain", Q::Normalized) = show(io, Q)


text:
using Test
using CompScienceMeshes

h = 2π / 51
Γ = meshcircle(1.0, h,2) ## creating the mesh first
points=vertices(Γ)     ## get the vertices
num_points=length(points)
segments=cellarray(Γ)  ## get the faces
Γ2 = barycentric_refinement(Γ) ## create the refinment of that one
points2=vertices(Γ2)     ## get the vertices
segments2=cellarray(Γ2)  ## get the faces

@test points[1] == points2[1]               # test if the first point in the original and refined mesh are the same
@test points[num_points]==points2[num_points]# test if the last point in the original has the same index in the refined one
@test segments[1,1] ==segments2[1,1]        # test if the face are starting with the same point in test
@test segments[1,2] ==segments2[2,2]        # test if the original second point in the first cell got separated by a single point

# test barycentric refinement of surface mesh
m = meshrectangle(1.0, 1.0, 0.25, 3)
f = barycentric_refinement(m)

@test numcells(f) == 6*numcells(m)

m1 = skeleton(m,1)
f1 = skeleton(f,1)
@test numcells(f1) == 2*numcells(m1) + 6*numcells(m)

m0 = skeleton(m,0)
@test numvertices(f) == numcells(m0) + numcells(m1) + numcells(m)

## test bisecting referinment of surfacic meshes
b = bisecting_refinement(m)
@test numcells(b) == 4*numcells(m)

b1 = skeleton(b,1)
@test numcells(b1) == 2*numcells(m1) + 3*numcells(m)
@test numvertices(b) == numcells(m0) + numcells(m1)

max_stars_repo_name: shakisparki/CompScienceMeshes

max_issues_repo_name: shakisparki/CompScienceMeshes

max_forks_repo_name: shakisparki/CompScienceMeshes


text:
using StructJuMP, MPI

numScens = 2
m = StructuredModel(numScens)

@defVar(m, x >= 1)
@defVar(m, y <= 2)

@addConstraint(m,  x + y == 1)
@addConstraint(m, -x + y <= 1)
setObjective(m, :Min, x*x + 0.5x*y + 0.25y*y - y)

rhs  = [5,4]
coef = [2,3]
qc   = [1,0.5]
ac   = [1,0.75]

@second_stage m scen begin
    bl = StructuredModel(parent=m)
    @defVar(bl, 0 <= w <= 1)
    @addConstraint(bl, coef[scen]w - x - y <= rhs[scen])
    @addConstraint(bl, coef[scen]w + x     == rhs[scen])
    setObjective(bl, :Min, qc[scen]*w*w + ac[scen]*w)
end


text:
# A mixed boundary-value problem described in the LaTeX file mixed.tex.

using FinElt
using FinElt.PlanarPoisson

h = 0.05
run(`gmsh -2 -clmax $h -o roundL.msh -format msh22 roundL.geo`)

mesh = read_msh_file("roundL.msh")
essential_bc = [ "Dirichlet" ]
vp = VariationalProblem(mesh, essential_bc)
add_bilin_form!(vp,  "Omega",     grad_dot_grad!)
add_lin_functnl!(vp, "Omega",     source_times_func!,     x->2.0)
add_lin_functnl!(vp, "NeumannSW", bdry_source_times_func!, x->-0.75)

A, b = assembled_linear_system(vp)

ufree = A \ b

u = complete_soln(ufree, vp)
write_pos_file("roundL.pos") do fid
    save_warp_nodal_scalar_field(u, "u", fid)
end

run(`gmsh mixed.script`)


text:
# FIXME: temporal patches to improve inferrability of Base methods
include("patches.jl")

# test utilities are extracted into a separate file for easier interactive testing from REPL
# i.e. julia-dev -i test/interactive_utils.jl
include("interactive_utils.jl")

# stuff used across tests
using Test

const FIXTURE_DIR = normpath(@__DIR__, "fixtures")

const ERROR_REPORTS_FROM_SUM_OVER_STRING = let
    interp, frame = analyze_call(sum, (String,))
    @test !isempty(interp.reports)
    interp.reports
end

function test_sum_over_string(ers)
    @test !isempty(ers)
    for target in ERROR_REPORTS_FROM_SUM_OVER_STRING
        @test any(ers) do er
            return er.msg == target.msg && er.sig == target.sig
        end
    end
end
test_sum_over_string(res::JET.VirtualProcessResult) =
    test_sum_over_string(res.inference_error_reports)
test_sum_over_string(interp::JETInterpreter) = test_sum_over_string(interp.reports)

return


text:
using MacroTools

# This is an incomplete implementation based on https://github.com/FluxML/Flux.jl/pull/299/files
# While it doesn't support all the features of numpy's einsum, it illustrates
# the important points about how we build functionality like this in Julia.

_permutedims(x, p) =
  p == collect(1:length(p)) ? x : # TODO 0.7
  p == [2, 1] ? :(transpose($x)) :
  :(permutedims($x, ($(p...),)))

_size(x, n) = Any[:(size($x, $i)) for i = 1:n]

_reshape(x, n, s) = # TODO use transpose
  s == _size(x, n) ? x :
  :(reshape($x, ($(s...),)))

function _expanddims(x, n, ds) # TODO use transpose
  is = _size(:x, n)
  foreach(d -> insert!(is, d, 1), ds)
  :(let x = $x
      reshape(x, ($(is...),))
    end)
end

function _squeezedims(x, n, ds)
  is = [d for (i, d) in enumerate(_size(:x, n)) if i ∉ ds]
  :(let x = $x
      reshape(x, ($(is...),))
    end)
end

function _einsum_pair(a, b, dims)
  (a, adims), (b, bdims) = a, b # TODO 0.7
  preserved = setdiff(intersect(adims, bdims), dims)
  broadcast = map(adims -> reduce(setdiff, (adims, preserved, dims)), (adims, bdims)) # TODO 0.7
  # TODO move preserved dims last
  aperm = sortperm(adims, by = i -> i in preserved ? -1 : i in broadcast[1] ? 0 : 1)
  bperm = sortperm(bdims, by = i -> i in preserved ? -1 : i in dims ? 0 : 1)
  a, b = _permutedims(a, aperm), _permutedims(b, bperm)
  adims, bdims = adims[aperm], bdims[bperm]
  if isempty(dims)
    b = _expanddims(b, length(bdims), length(preserved).+(1:length(broadcast[1])))
    :(broadcast(*, $a, $b)), vcat(adims[aperm], bdims[bperm][length(preserved)+1:end])
  else
    prod(xs) = isempty(xs) ? 0 : length(xs) == 1 ? xs[1] : :(prod(($(xs...),)))

    ashape = _size(a, length(adims))
    npreserve = prod(ashape[1:length(preserved)])
    aaxes = 1+length(preserved):length(adims)-length(dims)
    abroadcast = prod(ashape[aaxes])
    asum = prod(ashape[end-length(dims)+1:end])
    a = _reshape(a, length(adims), [abroadcast, asum]) # TODO preserve

    bshape = _size(b, length(bdims))
    bsum = prod(bshape[length(preserved)+1:end-length(broadcast[2])])
    baxes = length(bdims)-length(broadcast[2])+1:length(bdims)
    bbroadcast = prod(bshape[baxes])
    b = _reshape(b, length(bdims), [bsum, bbroadcast]) # TODO preserve

    ab = :($a*$b)
    shape = vcat(ashape[[1:length(preserved)..., aaxes...]], bshape[baxes])
    shape == [abroadcast, bbroadcast] || (ab = _reshape(ab, 2, shape))
    axes = vcat(adims[[1:length(preserved)..., aaxes...]], bdims[baxes])

    return ab, axes
  end
end

macro einsum(ex)
  @capture(ex, [out__] -> *(in__) | in_) || error("`@einsum [...] -> a[...] * b[...] * ...`")
  in isa Vector || (in = [in])
  # TODO rebinding, check dims
  in = map(in) do x
    @capture(x, a_[i__]) || error("Einsum input should be `a[i...]`, got `$x`")
    esc(a), i
  end
  all(length(unique(is)) == length(is) for (_, is) in in) || error("Diagonals not supported")
  labels = unique(vcat(map(x -> x[2], in)...))
  for i in labels
    count(in -> i ∈ in[2], in) > 2 && error("Not supported: index $i appears more than twice")
  end
  y = in[1]
  for i = 1:length(in)-1
    dims = setdiff(union(y[2], in[i+1][2]), out)
    y = _einsum_pair(y, in[i+1], dims)
  end
  reduce = setdiff(y[2], out)
  if !isempty(reduce)
    r = indexin(reduce, y[2])
    y = _squeezedims(:(sum($(y[1]), dims = ($(r...),))), length(y[2]), r),
        setdiff(y[2], reduce)
  end
  @assert sort(y[2]) == sort(out)
  return _permutedims(y[1], indexin(out, y[2]))
end

max_stars_repo_name: MikeInnes/zygote-paper

max_issues_repo_name: MikeInnes/zygote-paper

max_forks_repo_name: MikeInnes/zygote-paper