text:
using Distributed
using DistributedArrays
using DistributedArrays.SPMD
#@everywhere using SharedArrays
@everywhere using Distributed
@everywhere using DistributedArrays
@everywhere using DistributedArrays.SPMD

"""
The following function evaluates the electric field on a uniform grid from the electric potential.

    // Calculate electric field from potential
"""
function get_E_from_ϕ!(ϕ, E, dx)
      J = length(E)
      for j in 2:J-1
        E[j] = (ϕ[j-1] - ϕ[j+1]) / 2. / dx
      end
      E[1] = (ϕ[J] - ϕ[2]) / 2. / dx;
      E[J] = (ϕ[J-1] - ϕ[1]) / 2. / dx;
end

""" The following routine solves Poisson's equation in 1-D to find the instantaneous electric potential on a uniform grid.

// Solves 1-d Poisson equation:
//    d^u / dx^2 = v   for  0 <= x <= L
// Periodic boundary conditions:
//    u(x + L) = u(x),  v(x + L) = v(x)
// Arrays u and v assumed to be of length J.
// Now, jth grid point corresponds to
//    x_j = j dx  for j = 0,J-1
// where dx = L / J. L / (J-1) in Julia
// Also,
//    kappa = 2 pi / L
"""
function get_ϕ!(ϕ, ρ, κ)
  #V = fill(0.0+im*0.,J)
  #U = fill(0.0+im*0.,J÷2+1)
  J = length(ρ)
  # Fourier transform source term
  V = rfft(ρ)

  # Calculate Fourier transform of u

  V[1] =  0.
  for j in  2:(J÷2+1)
    V[j] = - V[j] / (j-1)^2 / κ^2
  end

  # Inverse Fourier transform to obtain u
  ϕ[:] = irfft(V,J)
end

"""
Takes out the mass of the grid function so that the sum is now null
"""
function filter_constant!(E)
  J = length(E)
  V = rfft(E)
  V[1] = 0.0 #extract the first component
  E[:] = irfft(V,J)
end

"""The routine below evaluates the electron number density on an evenly spaced mesh given the instantaneous electron coordinates.

// Evaluates electron number density n(0:J-1) from
// array r(0:N-1) of electron coordinates.
"""
function get_density_old_old!(r, n, dx)
  # Initialize
  N = length(r)
  J = length(n)
  fill!(n,0.0)
  # Evaluate number density.
  for i in 1:N
    j, y = get_index_and_distance(r[i],dx,L)
    if j < 1
        error("j = $j")
    end
    n[j] += (1. - y) / dx;
    if (j == J)
        n[1] += y / dx
    else
        n[j+1] += y / dx
    end
  end
end

function get_density_old!(u, n, p)
  L, N, J, κ, dx, order = p
  r = view(u,1:N)
  fill!(n,0.0)
  # Evaluate number density.
  for i in 1:N
    j, y = get_index_and_y(r[i],J,L)
    #if j < 1
    #    error("j = $j")
    #end
    if order == 1
      n[j] += (1. - y) / dx;
      if (j == J)
          n[1] += y / dx
      else
        n[j+1] += y / dx
      end
    end
    if order == 2
      if y <= 1/2
        n[j] += (3/4 - y^2) / dx
        if j == 1
          n[J] += (1-2y)^2/8 / dx
          n[2] += (1+2y)^2/8 / dx
        elseif j==J
          n[1] += (1+2y)^2/8 / dx
          n[J-1] += (1-2y)^2/8 / dx
        else
          n[j+1] += (1+2y)^2/8 / dx
          n[j-1] += (1-2y)^2/8 / dx
        end
      elseif y <= 3/2
        n[j] += (3 - 2y)^2/8 / dx
        if j == J-1
          n[J] += (3/4 - (1-y)^2) / dx
          n[1] += (1-2y)^2/8 / dx
        elseif j==J
          n[1] += (3/4 - (1-y)^2) / dx
          n[2] += (1-2y)^2/8 / dx
        else
          n[j+1] += (3/4 - (1-y)^2) / dx
          n[j+2] += (1-2y)^2/8 / dx
        end
      end
    end
  end
    #n .= n/n0 - 1.0 # return rho directly
end

function get_density!(u, n, p)
  L, N, J, κ, dx, order = p
  r = view(u,1:N)
  fill!(n,0.0)
  # Evaluate number density.
  for i in 1:N
    @inbounds j, y = get_index_and_y(r[i],J,L)
    for l in (-order):order
      @inbounds n[mod1(j + l, J)] += W(order, -y + l) / dx;
    end
  end
    #n .= n/n0 - 1.0 # return rho directly
end

function get_density_ro!(u, n, p)
  L, N, J, κ, dx, order = p
  r = view(u,1:2:2N-1)
  fill!(n,0.0)
  # Evaluate number density.
  for i in 1:N
    j, y = get_index_and_y(r[i],J,L)
    for l in (-order):order
      n[mod1(j + l, J)] += W(order, -y + l) / dx;
    end
  end
    #n .= n/n0 - 1.0 # return rho directly
end

"""The routine below evaluates the electron current on an evenly spaced mesh given the instantaneous electron coordinates.

// Evaluates electron number density S(0:J-1) from
// array r(0:N-1) of electron coordinates.
"""
function get_current_old!(u, S, p)
  L, N, J, κ, dx, order = p
  r = view(u,1:N)
  v = view(u,N+1:2N)
  fill!(S,0.0)
  # Evaluate number density.
  for i in 1:N
    #j, y = get_index_and_distance(r[i],dx,L)
    j, y = get_index_and_y(r[i],J,L)
    #if j < 1 || j > J
    #    error("j = $j")
    #end
    if order == 1
      S[j] += (1. - y)*v[i] / dx;
      if (j == J)
        S[1] += y*v[i] / dx
      else
        S[j+1] += y*v[i] / dx
      end
    end
    if order == 2
      if y <= 1/2
        S[j] += (3/4 - y^2) *v[i] / dx
        if j == 1
          S[J] += (1-2y)^2/8 *v[i] / dx
          S[2] += (1+2y)^2/8 *v[i] / dx
        elseif j==J
          S[1] += (1+2y)^2/8 *v[i] / dx
          S[J-1] += (1-2y)^2/8 *v[i] / dx
        else
          S[j+1] += (1+2y)^2/8 *v[i] / dx
          S[j-1] += (1-2y)^2/8 *v[i] / dx
        end
      elseif y <= 3/2
        S[j] += (3 - 2y)^2/8 *v[i] / dx
        if j == J-1
          S[J] += (3/4 - (1-y)^2) *v[i] / dx
          S[1] += (1-2y)^2/8 *v[i] / dx
        elseif j==J
          S[1] += (3/4 - (1-y)^2) *v[i] / dx
          S[2] += (1-2y)^2/8 *v[i] / dx
        else
          S[j+1] += (3/4 - (1-y)^2) *v[i] / dx
          S[j+2] += (1-2y)^2/8 *v[i] / dx
        end
      end
    end
  end
end

function get_current!(u, S, p)
  L, N, J, κ, dx, order = p
  r = view(u,1:N)
  v = view(u,N+1:2N)
  fill!(S,0.0)

  for i in 1:N
    @inbounds j, y = get_index_and_y(r[i],J,L)
    for l in (-order):order
      @inbounds S[mod1(j + l, J)] += W(order, -y + l) * v[i] / dx;
      #S[mod1(j + l, J)] += W_alt(order, -y + l) * v[i] / dx;
    end
  end
end

function get_current_ro!(u, S, p)
  L, N, J, κ, dx, order = p
  r = view(u,1:2:2N-1)
  v = view(u,2:2:2N)
  fill!(S,0.0)

  for i in 1:N
    @inbounds j, y = get_index_and_y(r[i],J,L)
    for l in (-order):order
      @inbounds S[mod1(j + l, J)] += W(order, -y + l) * v[i] / dx;
    end
  end
end

function get_current(u, S, p)
  L, N, J, κ, dx, order = p
  r = view(u,1:N)
  v = view(u,N+1:2N)
  fill!(S,0.0)

  for i in 1:N
    @inbounds j, y = get_index_and_y(r[i],J,L)
    for l in (-order):order
      @inbounds  S[mod1(j + l, J)] += W(order, -y + l) * v[i] / dx;
    end
  end
  return S[:]
end

function get_mini_current(pp, S, p)
  L, N, J, κ, dx, order = p
  fill!(S,0.0)
    j, y = get_index_and_y(pp.r,J,L)
    for l in (-order):order
      S[mod1(j + l, J)] = W(order, -y + l) * pp.v / dx;
    end
  return S[:]
end

function get_current_threads!(u, S, p)
  L, N, J, κ, dx, order, TS = p
  j = fill(Int64(1),nthreads())
  y = fill(Float64(1.0),nthreads())
  TS .= zeros(Float64)
  @threads for i in 1:N
    @inbounds j[threadid()], y[threadid()] = get_index_and_y(u[i], J, L)
    for l in (-order):-j[threadid()]
      @inbounds TS[J + j[threadid()] + l, threadid()] += W(order, -y[threadid()] + l) * u[N+i] / dx;
    end
    for l in max(-order,-j[threadid()]+1):min(order,J-j[threadid()])
      @inbounds TS[j[threadid()] + l, threadid()] += W(order, -y[threadid()] + l) * u[N+i] / dx;
    end
    for l in J-j[threadid()]+1:order
      @inbounds TS[j[threadid()] - J + l, threadid()] += W(order, -y[threadid()] + l) * u[N+i] / dx;
    end
    # for l in (-order):order
    #   @inbounds TS[mod1(j + l, J), threadid()] += W(order, -y + l) * v[i] / dx;
    # end
  end

  S .= zeros(Float64)
  @threads for i in 1:J
    for t in 1:nthreads()
      @inbounds S[i] += TS[i, t]
    end
  end
  S
end

function get_current_ro_par(Dpars::DArray, DS::DArray, p)
  L, N, J, κ, dx, order = p
  ran_pars = first(localindices(Dpars))
  #assuming ran_pars[1] is odd
  odd_range = first(ran_pars):2:last(ran_pars)
  LS = DS[:L]
  fill!(LS,0.0)
  for i in odd_range
    j, y = get_index_and_y(Dpars[i],J,L)
    for l in (-order):order
      LS[mod1(j + l, J),:] .+= W(order, -y + l) * Dpars[i+1] / dx;
    end
  end
end

"""
Calculates the RHS of the evolution equation.
Returns du
uses several functions which are passed as parameters
p = N, J, L, dx, Density!, Electric!, Poisson1D!
"""
function RHS(u,t,p_RHC)
    N, J, L, dx, n, E, ϕ, du, Density!, Electric!, Poisson1D! = p_RHC
    p = L, N, J, κ, dx
    get_density!(u[1:N], n, p)
    n0 = N/L
    get_ϕ!(ϕ, n/n0 .-1.0,κ)
    get_E_from_ϕ!(ϕ,E,dx)

    for i in 1:N
        #j, y = get_index_and_distance(u[i],dx,L)
        j, y = get_index_and_y(u[i],J,L)
        if (j == J)
            Efield = E[j] * (1. - y) + E[1] * y;
        else
            Efield = E[j] * (1. - y) + E[j+1] * y;
        end

        du[i] = u[N+i]
        du[N+i] = - Efield
    end
    return du[:]
end

"""
Calculates the RHS of the evolution equation.
Returns du
uses several functions which are passed as parameters
p = N, J, L, dx, n, S, du, get_current!
"""
function RHSC_old(u,t,p_RHSC)
  N, J, L, dx, order, n, S, du, get_density!, get_current! = p_RHSC
    p = L, N, J, κ, dx, order
    #get_density!(u, n, p)
    get_current!(u, S, p)
    E = view(u,2N+1:2N+J)
    n0 = N/L

    for i in 1:N
        #j, y = get_index_and_distance(u[i],dx,L)
        j, y = get_index_and_y(u[i],J,L)

        if order == 1
          if (j == J)
              Efield = E[j] * (1. - y) + E[1] * y;
          else
              Efield = E[j] * (1. - y) + E[j+1] * y;
          end
        elseif order == 2
          if y <= 1/2
            if j == 1
              Efield = E[j] * (3/4 - y^2) + E[2] * (1+2y)^2/8 + E[J] * (1-2y)^2/8
            elseif j==J
              Efield = E[j] * (3/4 - y^2) + E[1] * (1+2y)^2/8 + E[J-1] * (1-2y)^2/8
            else
              Efield = E[j] * (3/4 - y^2) + E[j+1] * (1+2y)^2/8 + E[j-1] * (1-2y)^2/8
            end
          elseif y <= 3/2
            if j == J-1
              Efield = E[j] * (3 - 2y)^2/8 + E[J] * (3/4 - (1-y)^2) + E[1] * (1-2y)^2/8
            elseif j == J
              Efield = E[j] * (3 - 2y)^2/8 + E[1] * (3/4 - (1-y)^2) + E[2] * (1-2y)^2/8
            else
              Efield = E[j] * (3 - 2y)^2/8 + E[j+1] * (3/4 - (1-y)^2) + E[j+2] * (1-2y)^2/8
            end
          end
        end

        du[i] = u[N+i]
        du[N+i] = - Efield
    end

      for j in 1:J
        du[2N+j] =  S[j]/n0 # particles have negative sign!
      end

    return du[:]
end

function RHSC(u,t,p_RHSC)
  if nthreads() == 1
    N, J, L, dx, order, n, S, du, get_density!, get_current!, Interpolate = p_RHSC
    p = L, N, J, κ, dx, order
    get_current!(u, S, p)
  else
    N, J, L, dx, order, n, S, du, get_density!, get_current!, Interpolate, TS = p_RHSC
    p = L, N, J, κ, dx, order, TS
    get_current_threads!(u, S, p)
  end

    #get_density!(u, n, p)

    E = view(u,2N+1:2N+J)
    n0 = N/L

    #@threads for i in 1:N
    for i in 1:N
        #j, y = get_index_and_y(u[i],J,L)
        #Efield = 0.0
        #for l in (-order):order
        #  Efield += E[mod1(j+l,J)] * W(order, -y + l)
        #end

        @inbounds du[i] = u[N+i]
        @inbounds du[N+i] = - Interpolate(order, E, u[i], J, L)
    end

    #@threads for j in 1:J
    for j in 1:J
        @inbounds du[2N+j] =  S[j]/n0 # particles have negative sign!
    end
    return du[:]
end

function Coordinate_test(r,L)
    if minimum(r) < 0.0
        error("negative coordinates")
    end
    if maximum(r) > L
        error("coordintates extend beyond L")
    end
end

function make_periodic_old!(r,L)
    N = length(r)
    for i in 1:N
      if (r[i] < 0.)
        r[i] = r[i] + L;
      end
      if (r[i] > L)
        r[i] = r[i] - L;
      end
    end
    return r[:]
end

function make_periodic!(r,L)
  return mod1.(r,L)
end

function RK4_Step!(f,y0,t0,h,p)
    k1 = h*f(y0,t0,p)
    k2 = h*f(y0+0.5*k1, t0+0.5*h,p)
    k3 = h*f(y0+0.5*k2, t0+0.5*h,p)
    k4 = h*f(y0+k3, t0+h,p)
    y0 .= y0 + (k1 + 2k2 + 2k3 + k4)/6
end

function get_index_and_distance(s,dx,L)
    #if s < 0
    #    s = s + L
    #end
    #if s > L
    #    s = s - L
    #end
    s = mod1(s,L)
    j = floor(Int64, s ÷ dx) + 1 #find the grid space where it is.
    if j > J || j < 1
      error("j = $j")
    end
    y = (s % dx)/dx #how far is there
    return j, y
end

"""
given a number s in between x_j and x_{j+1} computes y = (s - x_j)/dx and return j and y.
get_index_and_y(0.4,2,1)
dx = 0.5, j = 1, y = 0.4/0.5 = 0.8
get_index_and_y(0.7,2,1) = j = 2, 0.2/0.5= 0.4
"""
function get_index_and_y(s,J,L)
  s = (s/L*J + J)%J
  j = floor(Int,s) + 1
  #j = convert(Int64,s) + 1
  y = (s%1)
  return j, y
end

function get_energy(u,p)
  L, N, J = p
  dx = L/J
  n0 = N/L
  energy_K = 0.0
  energy_E = 0.0
  for i in 1:N
    energy_K = energy_K + u[N+i]^2
  end
  for j in 1:J
    energy_E = energy_E + u[2N+j]^2
  end

  return energy_K/2,  dx*energy_E /2 * n0
end

"""
Derivatives of shape functions
"""
function W(order::Int,y::Float64)
  y = abs(y)
  if order == 0
    return  (y <= 1/2) ? 1 : 0
  elseif order ==1
    return  (y <= 1) ? 1 - y : 0
  elseif order == 2
    return (y <= 1/2) ? 3/4 - y^2  : (((y > 1/2) && (y <= 3/2)) ? (3 - 2*y)^2 / 8 : 0)
  elseif order == 3
    return (y <= 1) ? 2/3 - y^2 + y^3 / 2 : (((y > 1) && (y <= 2)) ? (2 - y)^3 / 6 : 0)
  elseif order == 4
    return (y <= 1/2) ? 115/192 - 5y^2/8 + y^4/4 : (((y > 1/2) && (y <= 3/2)) ? (55 + 20y -120y^2 + 80y^3 - 16y^4)/96 : (((y > 3/2) && (y < 5/2)) ? (5 - 2y)^4/384 : 0))
  elseif order == 5
    return (y <= 1) ? 11/20 - y^2/2 + y^4/4 - y^5/12 : (((y > 1) && (y <= 2)) ? 17/40 + 5y/8 - 7y^2/4 + 5y^3/4 - 3y^4/8 + y^5/24 : (((y > 2) && (y < 3)) ? (3 - y)^5/120 : 0))
  else
    error("order = $order not yet implemented ")
  end
end

"""
Alternative definition
"""
function W_alt(order::Int,y::Float64)
  y = abs(y)
  if order == 0
    return  (y > 1/2) ? 0 : 1
  elseif order ==1
    return  (y > 1) ? 0 : 1 - y
  elseif order == 2
    #return (y <= 1/2) ? 3/4 - y^2  : (((y > 1/2) && (y <= 3/2)) ? (3 - 2*y)^2 / 8 : 0)
    return (y > 3/2) ? 0 : (((y > 1/2) && (y <= 3/2)) ? (3 - 2*y)^2 / 8 : 3/4 - y^2)
  elseif order == 3
    #return (y <= 1) ? 2/3 - y^2 + y^3 / 2 : (((y > 1) && (y <= 2)) ? (2 - y)^3 / 6 : 0)
    return (y > 2) ? 0 : (((y > 1) && (y <= 2)) ? (2 - y)^3 / 6 : 2/3 - y^2 + y^3 / 2)
  elseif order == 4
    #return (y <= 1/2) ? 115/192 - 5y^2/8 + y^4/4 : (((y > 1/2) && (y <= 3/2)) ? (55 + 20y -120y^2 + 80y^3 - 16y^4)/96 : (((y > 3/2) && (y < 5/2)) ? (5 - 2y)^4/384 : 0))
    return (y > 5/2) ? 0 : (((y > 3/2) && (y < 5/2)) ? (5 - 2y)^4/384 : (((y > 1/2) && (y <= 3/2)) ? (55 + 20y -120y^2 + 80y^3 - 16y^4)/96 : 115/192 - 5y^2/8 + y^4/4))
  elseif order == 5
    #return (y <= 1) ? 11/20 - y^2/2 + y^4/4 - y^5/12 : (((y > 1) && (y <= 2)) ? 17/40 + 5y/8 - 7y^2/4 + 5y^3/4 - 3y^4/8 + y^5/24 : (((y > 2) && (y < 3)) ? (3 - y)^5/120 : 0))
    return (y > 3) ? 0 : (((y > 2) && (y < 3)) ? (3 - y)^5/120 : (((y > 1) && (y <= 2)) ? 17/40 + 5y/8 - 7y^2/4 + 5y^3/4 - 3y^4/8 + y^5/24 : 11/20 - y^2/2 + y^4/4 - y^5/12))
  else
    error("order = $order not yet implemented ")
  end
end
"""
This are interpolation functions for getting the Electric field correct.
According the SHARP the second is better. Since it keeps momentum conservation.
"""
function Interpolate_1(order, vector, x, J, L)
  j, y = get_index_and_y(x,J,L)
  vi = 0.0
    for l in (-order+1):order
      vi += vector[mod1(j+l,J)] * W(order, -y + l)
    end
  return vi
end

function Interpolate_2(order, vector, x, J, L)
  j, y = get_index_and_y(x,J,L)
  vi = 0.0
    for l in (-order):order
      vi += (vector[mod1(j+l,J)] + vector[mod1(j+l+1,J)])/ 2 * W(order, -y + 1/2 + l)
    end
  return vi
end
"""
Terminado, pero sin probar.
"""
function Interpolate_per(order, vector, x, J, L)
  j, y = get_index_and_y(x,J,L)
  vi = 0.0
    for l in (-order):-j
      vi += vector[J+j+l] * W(order, -y + 1/2 + l)
    end
    for l in  max(-order,-j+1):min(order,J-j)
      vi += vector[j+l] * W(order, -y + 1/2 + l)
    end
    for l in J-j+1:order
      vi += vector[j+l] * W(order, -y + 1/2 + l)
    end
    # now we increase j by 1
    j += 1
    for l in (-order):-j
      vi += vector[J+j+l] * W(order, -y + 1/2 + l)
    end
    for l in  max(-order,-j+1):min(order,J-j)
      vi += vector[j+l] * W(order, -y + 1/2 + l)
    end
    for l in J-j+1:order
      vi += vector[j+l] * W(order, -y + 1/2 + l)
    end
  return vi / 2
end

"""
Structure and functions to work with particles
This is probably much faster for the position and
velocity are adjacent in memory.
"""

mutable struct Particles
  r ::Float64
  v ::Float64
end

#pars = Vector{Particles}(undef, N)
#par1 = par(0.0,0.0)

#pp = fill(par1,N)
"""
Makes particles out of a vector with positions first and then
velocities
"""

function make_particles!(par_dis, pars)
  N = length(par_dis)÷2
  for i in 1:N
      pars[i] = Particles(par_dis[i], par_dis[N+i])
  end
end

"""
reorder particle vector so that position and velocity are contiguous
"""
function reorder_particles!(u,uro)
  N = length(u)÷2
  for i in 1:N
    uro[2i - 1] = u[i]
    uro[2i] = u[N+i]
  end
end


text:
"""
    abstract type Triangulation{Dt,Dp}

A discredited physical domain associated with a `DiscreteModel{Dm,Dp}`.

`Dt` and `Dm` can be different.

The (mandatory) `Triangulation` interface can be tested with

- [`test_triangulation`](@ref)
"""
abstract type Triangulation{Dc,Dp} <: Grid{Dc,Dp} end

function get_background_model(t::Triangulation)
  @abstractmethod
end

function get_grid(t::Triangulation)
  @abstractmethod
end

# See possible types of glue below
function get_glue(t::Triangulation,::Val{d}) where d
  nothing
end

"""
    test_triangulation(trian::Triangulation)
"""
function test_triangulation(trian::Triangulation{Dc,Dp}) where {Dc,Dp}
  test_grid(trian)
  model = get_background_model(trian)
  @assert isa(model,DiscreteModel)
  grid = get_grid(trian)
  test_grid(grid)
end

# Grid interface
get_node_coordinates(trian::Triangulation) = get_node_coordinates(get_grid(trian))
get_cell_node_ids(trian::Triangulation) = get_cell_node_ids(get_grid(trian))
get_reffes(trian::Triangulation) = get_reffes(get_grid(trian))
get_cell_type(trian::Triangulation) = get_cell_type(get_grid(trian))
get_facet_normal(trian::Triangulation) = get_facet_normal(get_grid(trian))

# The following are not strictly needed, sine there is a default implementation for them.
# In any case, we delegate just in case the underlying grid defines more performant versions
get_cell_coordinates(trian::Triangulation) = get_cell_coordinates(get_grid(trian))
get_cell_ref_coordinates(trian::Triangulation) = get_cell_ref_coordinates(get_grid(trian))
get_cell_shapefuns(trian::Triangulation) = get_cell_shapefuns(get_grid(trian))
get_cell_map(trian::Triangulation) = get_cell_map(get_grid(trian))
get_cell_reffe(trian::Triangulation) = get_cell_reffe(get_grid(trian))
is_first_order(trian::Triangulation) = is_first_order(get_grid(trian))

# This is the most used glue, but others are possible, see e.g. SkeletonGlue.
struct FaceToFaceGlue{A,B,C}
  tface_to_mface::A
  tface_to_mface_map::B
  mface_to_tface::C
end

function is_change_possible(strian::Triangulation,ttrian::Triangulation)
  if strian === ttrian
    return true
  end
  @check get_background_model(strian) === get_background_model(ttrian) "Triangulations do not point to the same background discrete model!"
  D = num_cell_dims(strian)
  sglue = get_glue(strian,Val(D))
  tglue = get_glue(ttrian,Val(D))
  is_change_possible(sglue,tglue)
end

function is_change_possible(strian,ttrian)
  false
end

function is_change_possible(sglue::FaceToFaceGlue,tglue::FaceToFaceGlue)
  sglue.mface_to_tface != nothing
end

function best_target(trian1::Triangulation,trian2::Triangulation)
  @check is_change_possible(trian1,trian2)
  @check is_change_possible(trian2,trian1)
  D1 = num_cell_dims(trian1)
  D2 = num_cell_dims(trian2)
  glue1 = get_glue(trian1,Val(D2))
  glue2 = get_glue(trian2,Val(D1))
  best_target(trian1,trian2,glue1,glue2)
end

function best_target(
  trian1::Triangulation,trian2::Triangulation,glue1::FaceToFaceGlue,glue2::FaceToFaceGlue)
  # Return the background
  model = get_background_model(trian1)
  D = num_cell_dims(trian1)
  @assert num_cell_dims(trian2) == D
  Triangulation(ReferenceFE{D},model)
end

function get_active_model(t::Triangulation)
  compute_active_model(t)
end

function compute_active_model(t::Triangulation)
  D = num_cell_dims(t)
  glue = get_glue(t,Val(D))
  @assert glue.mface_to_tface !== nothing
  bgmodel = get_background_model(t)
  model = DiscreteModel(Polytope{D},bgmodel)
  _restrict_model(model,get_grid(t),glue.tface_to_mface)
end

function _restrict_model(model,grid::Grid,tface_to_mface)
  _restrict_model(model,tface_to_mface)
end

function _restrict_model(model,grid::GridPortion,tface_to_mface)
  @check grid.cell_to_parent_cell == tface_to_mface
  DiscreteModelPortion(model,grid)
end

function _restrict_model(model,tface_to_mface)
  DiscreteModelPortion(model,tface_to_mface)
end

function _restrict_model(model,tface_to_mface::IdentityVector)
  model
end

# This is the most basic Triangulation
# It represents a physical domain built using the faces of a DiscreteModel
struct BodyFittedTriangulation{Dt,Dp,A,B,C} <: Triangulation{Dt,Dp}
  model::A
  grid::B
  tface_to_mface::C
  function BodyFittedTriangulation(model::DiscreteModel,grid::Grid,tface_to_mface)
    Dp = num_point_dims(model)
    @assert Dp == num_point_dims(grid)
    Dt = num_cell_dims(grid)
    A = typeof(model)
    B = typeof(grid)
    C = typeof(tface_to_mface)
    new{Dt,Dp,A,B,C}(model,grid,tface_to_mface)
  end
end

get_background_model(trian::BodyFittedTriangulation) = trian.model
get_grid(trian::BodyFittedTriangulation) = trian.grid

function get_glue(trian::BodyFittedTriangulation{Dt},::Val{Dt}) where Dt
  tface_to_mface = trian.tface_to_mface
  tface_to_mface_map = Fill(GenericField(identity),num_cells(trian))
  if isa(tface_to_mface,IdentityVector) && num_faces(trian.model,Dt) == num_cells(trian)
    mface_to_tface = tface_to_mface
  else
    nmfaces = num_faces(trian.model,Dt)
    mface_to_tface = PosNegPartition(tface_to_mface,Int32(nmfaces))
  end
  FaceToFaceGlue(tface_to_mface,tface_to_mface_map,mface_to_tface)
end

#function get_glue(trian::BodyFittedTriangulation{Dt},::Val{Dm}) where {Dt,Dm}
#  @notimplemented
#end

function Base.view(t::BodyFittedTriangulation,ids::AbstractArray)
  model = t.model
  grid = view(t.grid,ids)
  tface_to_mface = lazy_map(Reindex(t.tface_to_mface),ids)
  BodyFittedTriangulation(model,grid,tface_to_mface)
end

get_triangulation(model) = Triangulation(model)

function Triangulation(
  ::Type{ReferenceFE{d}},model::DiscreteModel,filter::AbstractArray) where d
  mgrid = Grid(ReferenceFE{d},model)
  # Grid portion is OK here since this is usually used to
  # define a FE space
  tgrid = GridPortion(mgrid,filter)
  tface_to_mface = tgrid.cell_to_parent_cell
  BodyFittedTriangulation(model,tgrid,tface_to_mface)
end

function Triangulation(model::DiscreteModel,filter::AbstractArray)
  d = num_cell_dims(model)
  Triangulation(ReferenceFE{d},model,filter)
end

function Triangulation(
  ::Type{ReferenceFE{d}},
  model::DiscreteModel,
  labels::FaceLabeling;tags=nothing) where d

  if tags === nothing
    grid = Grid(ReferenceFE{d},model)
    tface_to_mface = IdentityVector(num_cells(grid))
    BodyFittedTriangulation(model,grid,tface_to_mface)
  else
    mface_to_mask = get_face_mask(labels,tags,d)
    Triangulation(ReferenceFE{d},model,mface_to_mask)
  end
end

function Triangulation(
  ::Type{ReferenceFE{d}},model::DiscreteModel;kwargs...) where d
  labels = get_face_labeling(model)
  Triangulation(ReferenceFE{d},model,labels;kwargs...)
end

function Triangulation(model::DiscreteModel;kwargs...)
  d = num_cell_dims(model)
  labels = get_face_labeling(model)
  Triangulation(ReferenceFE{d},model,labels;kwargs...)
end

function Triangulation(trian::Triangulation,args...;kwargs...)
  amodel = get_active_model(trian)
  dtrian = Triangulation(amodel,args...;kwargs...)
  CompositeTriangulation(trian,dtrian)
end

function Triangulation(trian::Triangulation)
  trian
end

function Triangulation(trian::Triangulation,x::AbstractArray{<:Integer})
  view(trian,x)
end

function Triangulation(trian::Triangulation,x::AbstractArray{<:Bool})
  y = findall(collect1d(x))
  view(trian,y)
end

function Interior(args...;kwargs...)
  Triangulation(args...;kwargs...)
end

# This is the low-level functionality to move from one Triangulation to another

function restrict(a::AbstractArray,b::AbstractArray)
  lazy_map(Reindex(a),b)
end

function extend(tface_to_val,mface_to_tface)
  @notimplemented
end

function extend(tface_to_val,mface_to_tface::IdentityVector)
  tface_to_val
end

function extend(tface_to_val,mface_to_tface::PosNegPartition)
  ipos_to_val, ineg_to_val = pos_neg_data(tface_to_val,mface_to_tface)
  i_to_iposneg = mface_to_tface
  lazy_map(PosNegReindex(ipos_to_val,ineg_to_val),i_to_iposneg)
end

function extend(a::LazyArray{<:Fill{typeof(transpose)}},b::PosNegPartition)
  c = a.args[1]
  d = extend(c,b)
  lazy_map(transpose,d)
end

function extend(a::LazyArray{<:Fill{typeof(linear_combination)}},b::PosNegPartition)
  d1 = extend(a.args[1],b)
  d2 = extend(a.args[2],b)
  lazy_map(linear_combination,d1,d2)
end

#function extend(a::LazyArray{<:Fill},b::PosNegPartition)
#  k = a.maps.value
#  args = map(i->extend(i,b),a.args)
#  lazy_map(k,args...)
#end

function pos_neg_data(
  ipos_to_val::AbstractArray{<:Number},i_to_iposneg::PosNegPartition)
  nineg = length(i_to_iposneg.ineg_to_i)
  ineg_to_val = Fill(zero(eltype(ipos_to_val)),nineg)
  ipos_to_val, ineg_to_val
end

function pos_neg_data(
  ipos_to_val::AbstractArray{<:AbstractArray{<:Number}},i_to_iposneg::PosNegPartition)
  nineg = length(i_to_iposneg.ineg_to_i)
  val = testitem(ipos_to_val)
  zs = 0 .* size(val)
  void = similar(val,eltype(val),zs)
  ineg_to_val = Fill(void,nineg)
  ipos_to_val, ineg_to_val
end

function pos_neg_data(
  ipos_to_val::AbstractArray{<:Field},i_to_iposneg::PosNegPartition)
  nineg = length(i_to_iposneg.ineg_to_i)
  ipos_to_v = lazy_map(VoidFieldMap(false),ipos_to_val)
  ineg_to_v = Fill(VoidField(testitem(ipos_to_val),true),nineg)
  ipos_to_v, ineg_to_v
end

function pos_neg_data(
  ipos_to_val::AbstractArray{<:AbstractArray{<:Field}},i_to_iposneg::PosNegPartition)
  _pos_neg_data_basis(ipos_to_val,i_to_iposneg)
end

function pos_neg_data(
  ipos_to_val::AbstractArray{<:AbstractArray{<:Dof}},i_to_iposneg::PosNegPartition)
  _pos_neg_data_basis(ipos_to_val,i_to_iposneg)
end

function _pos_neg_data_basis(ipos_to_val,i_to_iposneg)
  nineg = length(i_to_iposneg.ineg_to_i)
  ipos_to_v = lazy_map(VoidBasisMap(false),ipos_to_val)
  ineg_to_v = Fill(VoidBasis(testitem(ipos_to_val),true),nineg)
  ipos_to_v, ineg_to_v
end

function pos_neg_data(
  ipos_to_val::AbstractArray{<:ArrayBlock},i_to_iposneg::PosNegPartition)
  nineg = length(i_to_iposneg.ineg_to_i)
  val = testitem(ipos_to_val)
  void = _similar_empty(val)
  ineg_to_val = Fill(void,nineg)
  ipos_to_val, ineg_to_val
end

function pos_neg_data(
  ipos_to_val::AbstractArray{<:Tuple{<:Any,<:Any}},i_to_iposneg::PosNegPartition)
  nineg = length(i_to_iposneg.ineg_to_i)
  val = testitem(ipos_to_val)
  void = _similar_empty(val)
  ineg_to_val = Fill(void,nineg)
  ipos_to_val, ineg_to_val
end

function _similar_empty(val::AbstractArray)
  zs = 0 .* size(val)
  void = similar(val,eltype(val),zs)
end

function _similar_empty(val::ArrayBlock)
  a = deepcopy(val)
  for i in eachindex(a)
    if a.touched[i]
      a.array[i] = _similar_empty(a.array[i])
    end
  end
  a
end

function _similar_empty(val::Tuple)
  a, b = val
  a1 = _similar_empty(a)
  b1 = _similar_empty(b)
  (a1,b1)
end

# "Compose" triangulations

struct CompositeTriangulation{Dc,Dp,A,B} <: Triangulation{Dc,Dp}
  rtrian::A
  dtrian::B
  function CompositeTriangulation(
    rtrian::Triangulation, dtrian::Triangulation)
    @assert num_point_dims(rtrian) == num_point_dims(dtrian)
    Dp = num_point_dims(rtrian)
    Dc = num_cell_dims(dtrian)
    #@assert get_active_model(rtrian) === get_background_model(dtrian)
    A = typeof(rtrian)
    B = typeof(dtrian)
    new{Dc,Dp,A,B}(rtrian,dtrian)
  end
end

get_background_model(t::CompositeTriangulation) = get_background_model(t.rtrian)
get_active_model(t::CompositeTriangulation) = get_active_model(t.dtrian)
get_grid(t::CompositeTriangulation) = get_grid(t.dtrian)
get_facet_normal(t::CompositeTriangulation) = get_facet_normal(t.dtrian)
function get_glue(t::CompositeTriangulation,::Val{D}) where D
  Dr = num_cell_dims(t.rtrian)
  rglue = get_glue(t.rtrian,Val(D))
  dglue = get_glue(t.dtrian,Val(Dr))
  _compose_glues(rglue,dglue)
end

function _compose_glues(rglue,dglue)
  nothing
end

function _compose_glues(rglue::FaceToFaceGlue,dglue::FaceToFaceGlue)
  rface_to_mface = rglue.tface_to_mface
  dface_to_rface = dglue.tface_to_mface
  dface_to_mface = collect(lazy_map(Reindex(rface_to_mface),dface_to_rface))
  rface_to_mface_map = rglue.tface_to_mface_map
  dface_to_rface_map = dglue.tface_to_mface_map
  dface_to_mface_map1 = lazy_map(Reindex(rface_to_mface_map),dface_to_rface)
  dface_to_mface_map = lazy_map(∘,dface_to_mface_map1,dface_to_rface_map)
  mface_to_dface = nothing
  FaceToFaceGlue(dface_to_mface,dface_to_mface_map,mface_to_dface)
end

struct GenericTriangulation{Dc,Dp,A,B,C} <: Triangulation{Dc,Dp}
  grid::A
  model::B
  glue::C
  function GenericTriangulation(
    grid::Grid,
    model=nothing,
    glue=ntuple(i->nothing,num_cell_dims(grid)+1))
    Dc = num_cell_dims(grid)
    Dp = num_point_dims(grid)
    A = typeof(grid)
    B = typeof(model)
    C = typeof(glue)
    new{Dc,Dp,A,B,C}(grid,model,glue)
  end
end

get_grid(a::GenericTriangulation) = a.grid
get_glue(a::GenericTriangulation,::Val{D}) where D = a.glue[D+1]
function get_background_model(a::GenericTriangulation)
  @notimplementedif a.model === nothing "This triangulation object cannot be used to define a FE Space"
  a.model
end

struct TriangulationView{Dc,Dp,A,B} <: Triangulation{Dc,Dp}
  parent::A
  cell_to_parent_cell::B
  function TriangulationView(parent::Triangulation,cell_to_parent_cell::AbstractArray)
    Dc = num_cell_dims(parent)
    Dp = num_point_dims(parent)
    A = typeof(parent)
    B = typeof(cell_to_parent_cell)
    new{Dc,Dp,A,B}(parent,cell_to_parent_cell)
  end
end

Base.view(a::Triangulation,b::AbstractArray) = TriangulationView(a,b)

function TriangulationView(parent::Triangulation,parent_cell_to_mask::AbstractArray{Bool})
  cell_to_parent_cell = findall(collect1d(parent_cell_to_mask))
  TriangulationView(parent,cell_to_parent_cell)
end

function TriangulationView(parent::Triangulation,parent_cell_to_mask::AbstractVector{Bool})
  cell_to_parent_cell = findall(parent_cell_to_mask)
  TriangulationView(parent,cell_to_parent_cell)
end

function get_background_model(trian::TriangulationView)
  get_background_model(trian.parent)
end

function get_grid(trian::TriangulationView)
  view(get_grid(trian.parent),trian.cell_to_parent_cell)
end

function get_glue(t::TriangulationView,::Val{d}) where d
  parent = get_glue(t.parent,Val(d))
  parent === nothing && return nothing
  view(parent,t.cell_to_parent_cell)
end

function Base.view(glue::FaceToFaceGlue,ids::AbstractArray)
  tface_to_mface = lazy_map(Reindex(glue.tface_to_mface),ids)
  tface_to_mface_map = lazy_map(Reindex(glue.tface_to_mface_map),ids)
  if glue.mface_to_tface === nothing
    mface_to_tface = nothing
  else
    nmfaces = length(glue.mface_to_tface)
    mface_to_tface = PosNegPartition(tface_to_mface,Int32(nmfaces))
  end
  FaceToFaceGlue(tface_to_mface,tface_to_mface_map,mface_to_tface)
end

function get_facet_normal(trian::TriangulationView)
  n = get_facet_normal(trian.parent)
  restrict(n,trian.cell_to_parent_cell)
end

function get_cell_map(trian::TriangulationView)
  lazy_map(Reindex(get_cell_map(trian.parent)),trian.cell_to_parent_cell)
end

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    NamedTuple

`NamedTuple`s are, as their name suggests, named [`Tuple`](@ref)s. That is, they're a
tuple-like collection of values, where each entry has a unique name, represented as a
[`Symbol`](@ref). Like `Tuple`s, `NamedTuple`s are immutable; neither the names nor the values
can be modified in place after construction.

Accessing the value associated with a name in a named tuple can be done using field
access syntax, e.g. `x.a`, or using [`getindex`](@ref), e.g. `x[:a]`. A tuple of the
names can be obtained using [`keys`](@ref), and a tuple of the values can be obtained
using [`values`](@ref).

!!! note
    Iteration over `NamedTuple`s produces the *values* without the names. (See example
    below.) To iterate over the name-value pairs, use the [`pairs`](@ref) function.

The [`@NamedTuple`](@ref) macro can be used for conveniently declaring `NamedTuple` types.

# Examples
```jldoctest
julia> x = (a=1, b=2)
(a = 1, b = 2)

julia> x.a
1

julia> x[:a]
1

julia> keys(x)
(:a, :b)

julia> values(x)
(1, 2)

julia> collect(x)
2-element Array{Int64,1}:
 1
 2

julia> collect(pairs(x))
2-element Array{Pair{Symbol,Int64},1}:
 :a => 1
 :b => 2
```

In a similar fashion as to how one can define keyword arguments programmatically,
a named tuple can be created by giving a pair `name::Symbol => value` or splatting
an iterator yielding such pairs after a semicolon inside a tuple literal:

```jldoctest
julia> (; :a => 1)
(a = 1,)

julia> keys = (:a, :b, :c); values = (1, 2, 3);

julia> (; zip(keys, values)...)
(a = 1, b = 2, c = 3)
```
"""
Core.NamedTuple

if nameof(@__MODULE__) === :Base

@eval function NamedTuple{names,T}(args::Tuple) where {names, T <: Tuple}
    if length(args) != length(names::Tuple)
        throw(ArgumentError("Wrong number of arguments to named tuple constructor."))
    end
    # Note T(args) might not return something of type T; e.g.
    # Tuple{Type{Float64}}((Float64,)) returns a Tuple{DataType}
    $(Expr(:splatnew, :(NamedTuple{names,T}), :(T(args))))
end

function NamedTuple{names}(nt::NamedTuple) where {names}
    if @generated
        types = Tuple{(fieldtype(nt, n) for n in names)...}
        Expr(:new, :(NamedTuple{names, $types}), Any[ :(getfield(nt, $(QuoteNode(n)))) for n in names ]...)
    else
        types = Tuple{(fieldtype(typeof(nt), n) for n in names)...}
        NamedTuple{names, types}(Tuple(getfield(nt, n) for n in names))
    end
end

NamedTuple{names, T}(itr) where {names, T <: Tuple} = NamedTuple{names, T}(T(itr))
NamedTuple{names}(itr) where {names} = NamedTuple{names}(Tuple(itr))

end # if Base

length(t::NamedTuple) = nfields(t)
iterate(t::NamedTuple, iter=1) = iter > nfields(t) ? nothing : (getfield(t, iter), iter + 1)
firstindex(t::NamedTuple) = 1
lastindex(t::NamedTuple) = nfields(t)
getindex(t::NamedTuple, i::Int) = getfield(t, i)
getindex(t::NamedTuple, i::Symbol) = getfield(t, i)
indexed_iterate(t::NamedTuple, i::Int, state=1) = (getfield(t, i), i+1)
isempty(::NamedTuple{()}) = true
isempty(::NamedTuple) = false
empty(::NamedTuple) = NamedTuple()

convert(::Type{NamedTuple{names,T}}, nt::NamedTuple{names,T}) where {names,T<:Tuple} = nt
convert(::Type{NamedTuple{names}}, nt::NamedTuple{names}) where {names} = nt

function convert(::Type{NamedTuple{names,T}}, nt::NamedTuple{names}) where {names,T<:Tuple}
    NamedTuple{names,T}(T(nt))
end

if nameof(@__MODULE__) === :Base
    Tuple(nt::NamedTuple) = (nt...,)
    (::Type{T})(nt::NamedTuple) where {T <: Tuple} = convert(T, Tuple(nt))
end

function show(io::IO, t::NamedTuple)
    n = nfields(t)
    for i = 1:n
        # if field types aren't concrete, show full type
        if typeof(getfield(t, i)) !== fieldtype(typeof(t), i)
            show(io, typeof(t))
            print(io, "(")
            show(io, Tuple(t))
            print(io, ")")
            return
        end
    end
    if n == 0
        print(io, "NamedTuple()")
    else
        typeinfo = get(io, :typeinfo, Any)
        print(io, "(")
        for i = 1:n
            print(io, fieldname(typeof(t),i), " = ")
            show(IOContext(io, :typeinfo =>
                           t isa typeinfo <: NamedTuple ? fieldtype(typeinfo, i) : Any),
                 getfield(t, i))
            if n == 1
                print(io, ",")
            elseif i < n
                print(io, ", ")
            end
        end
        print(io, ")")
    end
end

eltype(::Type{NamedTuple{names,T}}) where {names,T} = eltype(T)

==(a::NamedTuple{n}, b::NamedTuple{n}) where {n} = Tuple(a) == Tuple(b)
==(a::NamedTuple, b::NamedTuple) = false

isequal(a::NamedTuple{n}, b::NamedTuple{n}) where {n} = isequal(Tuple(a), Tuple(b))
isequal(a::NamedTuple, b::NamedTuple) = false

_nt_names(::NamedTuple{names}) where {names} = names
_nt_names(::Type{T}) where {names,T<:NamedTuple{names}} = names

hash(x::NamedTuple, h::UInt) = xor(objectid(_nt_names(x)), hash(Tuple(x), h))

isless(a::NamedTuple{n}, b::NamedTuple{n}) where {n} = isless(Tuple(a), Tuple(b))
# TODO: case where one argument's names are a prefix of the other's

same_names(::NamedTuple{names}...) where {names} = true
same_names(::NamedTuple...) = false

# NOTE: this method signature makes sure we don't define map(f)
function map(f, nt::NamedTuple{names}, nts::NamedTuple...) where names
    if !same_names(nt, nts...)
        throw(ArgumentError("Named tuple names do not match."))
    end
    NamedTuple{names}(map(f, map(Tuple, (nt, nts...))...))
end

# a version of `in` for the older world these generated functions run in
@pure function sym_in(x::Symbol, itr::Tuple{Vararg{Symbol}})
    for y in itr
        y === x && return true
    end
    return false
end

@pure function merge_names(an::Tuple{Vararg{Symbol}}, bn::Tuple{Vararg{Symbol}})
    names = Symbol[an...]
    for n in bn
        if !sym_in(n, an)
            push!(names, n)
        end
    end
    (names...,)
end

@pure function merge_types(names::Tuple{Vararg{Symbol}}, a::Type{<:NamedTuple}, b::Type{<:NamedTuple})
    bn = _nt_names(b)
    Tuple{Any[ fieldtype(sym_in(n, bn) ? b : a, n) for n in names ]...}
end

"""
    merge(a::NamedTuple, bs::NamedTuple...)

Construct a new named tuple by merging two or more existing ones, in a left-associative
manner. Merging proceeds left-to-right, between pairs of named tuples, and so the order of fields
present in both the leftmost and rightmost named tuples take the same position as they are found in the
leftmost named tuple. However, values are taken from matching fields in the rightmost named tuple that
contains that field. Fields present in only the rightmost named tuple of a pair are appended at the end.
A fallback is implemented for when only a single named tuple is supplied,
with signature `merge(a::NamedTuple)`.

!!! compat "Julia 1.1"
    Merging 3 or more `NamedTuple` requires at least Julia 1.1.

# Examples
```jldoctest
julia> merge((a=1, b=2, c=3), (b=4, d=5))
(a = 1, b = 4, c = 3, d = 5)
```

```jldoctest
julia> merge((a=1, b=2), (b=3, c=(d=1,)), (c=(d=2,),))
(a = 1, b = 3, c = (d = 2,))
```
"""
function merge(a::NamedTuple{an}, b::NamedTuple{bn}) where {an, bn}
    if @generated
        names = merge_names(an, bn)
        types = merge_types(names, a, b)
        vals = Any[ :(getfield($(sym_in(n, bn) ? :b : :a), $(QuoteNode(n)))) for n in names ]
        :( NamedTuple{$names,$types}(($(vals...),)) )
    else
        names = merge_names(an, bn)
        types = merge_types(names, typeof(a), typeof(b))
        NamedTuple{names,types}(map(n->getfield(sym_in(n, bn) ? b : a, n), names))
    end
end

merge(a::NamedTuple{()}, b::NamedTuple) = b

merge(a::NamedTuple, b::Iterators.Pairs{<:Any,<:Any,<:Any,<:NamedTuple}) = merge(a, b.data)

merge(a::NamedTuple, b::NamedTuple, cs::NamedTuple...) = merge(merge(a, b), cs...)

merge(a::NamedTuple) = a

"""
    merge(a::NamedTuple, iterable)

Interpret an iterable of key-value pairs as a named tuple, and perform a merge.

```jldoctest
julia> merge((a=1, b=2, c=3), [:b=>4, :d=>5])
(a = 1, b = 4, c = 3, d = 5)
```
"""
function merge(a::NamedTuple, itr)
    names = Symbol[]
    vals = Any[]
    inds = IdDict{Symbol,Int}()
    for (k::Symbol, v) in itr
        oldind = get(inds, k, 0)
        if oldind > 0
            vals[oldind] = v
        else
            push!(names, k)
            push!(vals, v)
            inds[k] = length(names)
        end
    end
    merge(a, NamedTuple{(names...,)}((vals...,)))
end

keys(nt::NamedTuple{names}) where {names} = names
values(nt::NamedTuple) = Tuple(nt)
haskey(nt::NamedTuple, key::Union{Integer, Symbol}) = isdefined(nt, key)
get(nt::NamedTuple, key::Union{Integer, Symbol}, default) = haskey(nt, key) ? getfield(nt, key) : default
get(f::Callable, nt::NamedTuple, key::Union{Integer, Symbol}) = haskey(nt, key) ? getfield(nt, key) : f()
tail(t::NamedTuple{names}) where names = NamedTuple{tail(names)}(t)

@pure function diff_names(an::Tuple{Vararg{Symbol}}, bn::Tuple{Vararg{Symbol}})
    names = Symbol[]
    for n in an
        if !sym_in(n, bn)
            push!(names, n)
        end
    end
    (names...,)
end

"""
    structdiff(a::NamedTuple{an}, b::Union{NamedTuple{bn},Type{NamedTuple{bn}}}) where {an,bn}

Construct a copy of named tuple `a`, except with fields that exist in `b` removed.
`b` can be a named tuple, or a type of the form `NamedTuple{field_names}`.
"""
function structdiff(a::NamedTuple{an}, b::Union{NamedTuple{bn}, Type{NamedTuple{bn}}}) where {an, bn}
    if @generated
        names = diff_names(an, bn)
        types = Tuple{Any[ fieldtype(a, n) for n in names ]...}
        vals = Any[ :(getfield(a, $(QuoteNode(n)))) for n in names ]
        :( NamedTuple{$names,$types}(($(vals...),)) )
    else
        names = diff_names(an, bn)
        types = Tuple{Any[ fieldtype(typeof(a), n) for n in names ]...}
        NamedTuple{names,types}(map(n->getfield(a, n), names))
    end
end

"""
    setindex(nt::NamedTuple, val, key::Symbol)

Constructs a new `NamedTuple` with the key `key` set to `val`.
If `key` is already in the keys of `nt`, `val` replaces the old value.

```jldoctest
julia> nt = (a = 3,)
(a = 3,)

julia> Base.setindex(nt, 33, :b)
(a = 3, b = 33)

julia> Base.setindex(nt, 4, :a)
(a = 4,)

julia> Base.setindex(nt, "a", :a)
(a = "a",)
```
"""
function setindex(nt::NamedTuple, v, idx::Symbol)
    merge(nt, (; idx => v))
end

"""
    @NamedTuple{key1::Type1, key2::Type2, ...}
    @NamedTuple begin key1::Type1; key2::Type2; ...; end

This macro gives a more convenient syntax for declaring `NamedTuple` types. It returns a `NamedTuple`
type with the given keys and types, equivalent to `NamedTuple{(:key1, :key2, ...), Tuple{Type1,Type2,...}}`.
If the `::Type` declaration is omitted, it is taken to be `Any`.   The `begin ... end` form allows the
declarations to be split across multiple lines (similar to a `struct` declaration), but is otherwise
equivalent.

For example, the tuple `(a=3.1, b="hello")` has a type `NamedTuple{(:a, :b),Tuple{Float64,String}}`, which
can also be declared via `@NamedTuple` as:

```jldoctest
julia> @NamedTuple{a::Float64, b::String}
NamedTuple{(:a, :b),Tuple{Float64,String}}

julia> @NamedTuple begin
           a::Float64
           b::String
       end
NamedTuple{(:a, :b),Tuple{Float64,String}}
```

!!! compat "Julia 1.5"
    This macro is available as of Julia 1.5.
"""
macro NamedTuple(ex)
    Meta.isexpr(ex, :braces) || Meta.isexpr(ex, :block) ||
        throw(ArgumentError("@NamedTuple expects {...} or begin...end"))
    decls = filter(e -> !(e isa LineNumberNode), ex.args)
    all(e -> e isa Symbol || Meta.isexpr(e, :(::)), decls) ||
        throw(ArgumentError("@NamedTuple must contain a sequence of name or name::type expressions"))
    vars = [QuoteNode(e isa Symbol ? e : e.args[1]) for e in decls]
    types = [esc(e isa Symbol ? :Any : e.args[2]) for e in decls]
    return :(NamedTuple{($(vars...),), Tuple{$(types...)}})
end

max_stars_repo_name: fabmazz/julia

max_issues_repo_name: fabmazz/julia

max_forks_repo_name: fabmazz/julia


text:
export ArrayComparer, fn, reset!

type ArrayComparer <: AbstractComparer
    c::Comparer
end

ArrayComparer() = ArrayComparer(Comparer())

ArrayComparer(fn::Function) = ArrayComparer(Comparer(fn))

fn(ac::ArrayComparer) = (a::Vector, b::Vector) ->
begin
    baseFn = fn(ac.c)

    all = [map(i -> (i, 1), a); map(i -> (i, 2), b)]

    sort!(all, by = i -> i[1], lt = baseFn)

    score = 0
    arr1Points = 0
    arr2Points = 0

    for idx in 1:length(all)
        if idx > 1 && baseFn(all[idx - 1][1], all[idx][1])
            score += 1
        end
        if all[idx][2] == 1
            arr1Points += score
        else
            arr2Points += score
        end
    end

    arr1Points < arr2Points
end

function reset!(ac::ArrayComparer)
    reset!(ac.c)
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
const DOC_FINITE =
    "`AbstractArray{<:Finite}` (multiclass classification)"
const DOC_FINITE_BINARY =
    "`AbstractArray{<:Finite{2}}` (binary classification)"
const DOC_ORDERED_FACTOR =
    "`AbstractArray{<:OrderedFactor}` (classification of ordered target)"
const DOC_ORDERED_FACTOR_BINARY =
    "`AbstractArray{<:OrderedFactor{2}}` "*
    "(binary classification where choice of \"true\" effects the measure)"
const DOC_CONTINUOUS = "`AbstractArray{Continuous}` (regression)"

## TRAITS

is_measure_type(::Any) = false

# The following traits, with fallbacks defined in
# StatisticalTraits.jl, make sense for some or all measures:

const MEASURE_TRAITS = [:name,
                        :instances,
                        :human_name,
                        :target_scitype,
                        :supports_weights,
                        :supports_class_weights,
                        :prediction_type,
                        :orientation,
                        :reports_each_observation,
                        :aggregation,
                        :is_feature_dependent,
                        :docstring,
                        :distribution_type]

## FOR BUILT-IN MEASURES (subtyping Measure)

abstract type Measure <: MLJType end
is_measure_type(::Type{<:Measure}) = true
is_measure(m) = is_measure_type(typeof(m))

# docstring fall-back:
_decorate(s::AbstractString) = "`$s`"
_decorate(v::Vector{<:AbstractString}) = join(_decorate.(v), ", ")
function MMI.docstring(M::Type{<:Measure})
    list = _decorate(instances(M))
    ret = "`$(name(M))` - $(human_name(M)) type"
    isempty(list) || (ret *= " with instances $list")
    ret *= ". "
    return ret
end

# display:
show_as_constructed(::Type{<:Measure}) = true

# info (see also src/init.jl):
function MLJScientificTypes.info(M, ::Val{:measure_type})
    values = Tuple(@eval($trait($M)) for trait in MEASURE_TRAITS)
    return NamedTuple{Tuple(MEASURE_TRAITS)}(values)
end
MLJScientificTypes.info(m, ::Val{:measure}) = info(typeof(m))

## AGGREGATION

(::Sum)(v) = sum(v)
(::Sum)(v::LittleDict) = sum(values(v))

(::Mean)(v) = mean(v)
(::Mean)(v::LittleDict) = mean(values(v))

(::RootMeanSquare)(v) = sqrt(mean(v.^2))

aggregate(v, measure) = aggregation(measure)(v)

# aggregation is no-op on scalars:
const MeasureValue = Union{Real,Tuple{<:Real,<:Real}} # number or interval
aggregate(x::MeasureValue, measure) = x

## DISPATCH FOR EVALUATION

# yhat - predictions (point or probabilisitic)
# X - features
# y - target observations
# w - per-observation weights

function value(measure, yhat, X, y, w)
    vfdep     = Val(is_feature_dependent(measure))
    vsweights = Val(supports_weights(measure) ||
                    supports_class_weights(measure))
    return value(measure, yhat, X, y, w, vfdep, vsweights)
end

## DEFAULT EVALUATION INTERFACE

#  is feature independent, weights not supported:
value(m, yhat, X, y, w, ::Val{false}, ::Val{false}) = m(yhat, y)

#  is feature dependent:, weights not supported:
value(m, yhat, X, y, w, ::Val{true}, ::Val{false}) = m(yhat, X, y)

#  is feature independent, weights supported:
value(m, yhat, X, y, w,         ::Val{false}, ::Val{true}) = m(yhat, y, w)
value(m, yhat, X, y, ::Nothing, ::Val{false}, ::Val{true}) = m(yhat, y)

#  is feature dependent, weights supported:
value(m, yhat, X, y, w,         ::Val{true}, ::Val{true}) = m(yhat, X, y, w)
value(m, yhat, X, y, ::Nothing, ::Val{true}, ::Val{true}) = m(yhat, X, y)

## helper

function check_pools(ŷ, y)
    levels(y) == levels(ŷ[1]) ||
        error("Conflicting categorical pools found "*
              "in observations and predictions. ")
    return nothing
end

## INCLUDE SPECIFIC MEASURES AND TOOLS

include("meta_utilities.jl")
include("continuous.jl")
include("confusion_matrix.jl")
include("finite.jl")
include("loss_functions_interface.jl")

## DEFAULT MEASURES
default_measure(T, S) = nothing

# Deterministic + Continuous / Count ==> RMS
default_measure(::Type{<:Deterministic},
                ::Type{<:Union{Vec{<:Continuous}, Vec{<:Count}}}) = rms

# Deterministic + Finite ==> Misclassification rate
default_measure(::Type{<:Deterministic},
                ::Type{<:Vec{<:Finite}}) = misclassification_rate

# default_measure(::Type{Probabilistic},
#                 ::Type{<:Union{Vec{<:Continuous},
#                                Vec{<:Count}}}) = ???

# Probabilistic + Finite ==> Cross entropy
default_measure(::Type{<:Probabilistic},
                ::Type{<:Vec{<:Finite}}) = cross_entropy

# Fallbacks
default_measure(M::Type{<:Supervised}) = default_measure(M, target_scitype(M))
default_measure(::M) where M <: Supervised = default_measure(M)


text:
export tauchen_2d!,
    tauchen_1d!,
    tauchen_2d,
    tauchen_1d,
    minp_default

function bvn_upperlower_cdf(xlim, ylim, r)
    xl,xu = xlim
    yl,yu = ylim
    xl < xu && yl < yu || throw(DomainError())
    return bvncdf(xu,yu,r) - bvncdf(xl,yu,r) - bvncdf(xu,yl,r) + bvncdf(xl,yl,r)
end

function plus_minus_dx(x::Number, xspace::AbstractRange, scaling::Number=1)
  dx = step(xspace)/2
  x == first(xspace) && return (-Inf, dx,).*scaling
  x == last(xspace)  && return (-dx, Inf,).*scaling
  return (-dx, dx).*scaling
end

x_pm_Δ(x::Number, xspace::AbstractVector) = x .+ plus_minus_dx(x,xspace)

function tauchen_1d!(P::AbstractMatrix, S::AbstractVector, μ::Function, σ2::Number)
    σ2 > 0 || throw(DomainError())
    all(length(S) .== size(P)) || throw(DimensionMismatch())

    sigma = sqrt(σ2)
    invsigma = inv(sigma)

    for (j,s1) in enumerate(S)
      pm_dz = plus_minus_dx(s1, S, invsigma)
      for (i,s0) in enumerate(S)
        zm, zp = ((s1 - μ(s0))/sigma) .+ pm_dz
        P[i,j] = normcdf(zp) - normcdf(zm)
      end
    end
end

function tauchen_1d(S::AbstractVector{T}, μ::Function, σ2::Number) where {T<:Real}
    n = length(S)
    P = Matrix{T}(undef,n,n)
    tauchen_1d!(P,S,μ,σ2)
    return P
end

function tauchen_2d(S::Base.Iterators.ProductIterator, μ::Function, Σ::AbstractMatrix{T}) where {T<:Real}
  n = length(S)
  P = Matrix{T}(undef,n,n)
  tauchen_2d!(P,S,μ,Σ)
  return P
end

function tauchen_2d!(P::AbstractMatrix, S::Base.Iterators.ProductIterator, μ::Function, Σ::AbstractMatrix)
    ndims(S) == 2
    2 == checksquare(Σ) || throw(DimensionMismatch("Σ must be 2x2"))
    issymmetric(Σ) || throw(error("Σ must be symmetric"))
    length(S) == checksquare(P) || throw(DimensionMismatch())

    sigmas = sqrt.(diag(Σ))
    invsigmas = inv.(sigmas)
    ρ = Σ[1,2] / prod(sigmas)

    for (j,s1) in enumerate(S)
      pm_dz = plus_minus_dx.(s1, S.iterators, invsigmas)
      for (i,s0) in enumerate(S)
        z = (s1 .- μ(s0))./sigmas
        P[i,j] = bvn_upperlower_cdf(z[1] .+ pm_dz[1], z[2] .+ pm_dz[2], ρ)
      end
    end
end

minp_default() = 1e-5

function zero_out_small_probs!(Q::AbstractMatrix, P::AbstractMatrix, minp::Real)
    checksquare(P) == checksquare(Q) || throw(DimensionMismatch())
    @assert all(sum(P, dims=2) .≈ 1)
    for i in eachindex(P)
        p = P[i]
        Q[i] = p > minp ? p : 0
    end
    Q ./= sum(Q, dims=2)
    return Q
end

zero_out_small_probs!(P,minp) = zero_out_small_probs!(P,P,minp)
zero_out_small_probs( P,minp) = zero_out_small_probs!(similar(P), P, minp)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SPIRV_LLVM_Translator_unified_jll
using Base
using Base: UUID
using LazyArtifacts
Base.include(@__MODULE__, joinpath("..", ".pkg", "platform_augmentation.jl"))
import JLLWrappers

JLLWrappers.@generate_main_file_header("SPIRV_LLVM_Translator_unified")
JLLWrappers.@generate_main_file("SPIRV_LLVM_Translator_unified", UUID("85f0d8ed-5b39-5caa-b1ae-7472de402361"))
end  # module SPIRV_LLVM_Translator_unified_jll


text:
function weak_phs_FEM(Ne,N1, N2,a,b, weightQ1, weightQ2)
	# this function uses the partioned FEM method, finding a single PHS model
	# for the whole (1D) domain.
	# TO DO: - correct the input matrix, and allow for higher-order polynomials
	#	N1  (old = Norder+1)  variables x1, e1 and v1
	#	N2  (old = Norder)  #    variables x2, e2 and v2

		x1,w1,P = lglnodes(N1-1,(b-a)/Ne,0)    # discretization of x1 variables
		x1i,w1i,Pi = lglnodes(N1,0,(b-a)/Ne)  # these points are used to integrate the mass matrix using quadrature
		if N2 > 1
			x2,w2 = lglnodes(N2-1,(b-a)/Ne,0)
		else
			x2,w2 = [(a+b)/2], [(b-a)/Ne];
		end
		#M = massmatrix(ei,xi,we);
		#xquad, wquad = lgwt(Ne,a,b)\
		#M = massmatrix(ei,xi, xquad, wquad);
		M1 = massmatrix(x1,x1, x1i,w1i);
		x1i2,w1i2,Pi2 = lglnodes(N1*4,0,(b-a)/Ne)  # these points are used to integrate the mass matrix using quadrature
		L1 = ELid -> massmatrix(x1,x1, x1i2,w1i2, (x, zi, i) -> PortHamiltonian.leg_pol(x,zi,i)*weightQ1(x+(ELid-1)*(b-a)/Ne));

#		if N2 == 1
#			M2 = massmatrix(x2,x2, w2);
#		else
		x2i,w2i,P2i = lglnodes(N2,0,(b-a)/Ne)  # these points are used to integrate the mass matrix using quadrature
		M2 = massmatrix(x2,x2, x2i,w2i);
		x2i2,w2i2,P2i2 = lglnodes(N2*4,0,(b-a)/Ne)
		L2 = ELid -> massmatrix(x2,x2, x2i2,w2i2, (x, zi, i) -> PortHamiltonian.leg_pol(x,zi,i)*weightQ2(x+(ELid-1)*(b-a)/Ne));
#		end

		D = dermatrix(x1,x2,1)'*M2 ;

		p0 = map(i->leg_pol(0,x1,i), 1:length(x1))
		pL = map(i->leg_pol((b-a)/Ne,x1,i), 1:length(x1))
		B = [p0 pL]

		M1full = zeros(Ne*(N1-1)+1,Ne*(N1-1)+1);
		L1full = zeros(Ne*(N1-1)+1,Ne*(N1-1)+1);
		if N2 == 1
			M2full = zeros(Ne,Ne);
			L2full = zeros(Ne,Ne);
		else
			M2full = zeros((Ne)*(N2-1)+1,(Ne)*(N2-1)+1);
			L2full = zeros((Ne)*(N2-1)+1,(Ne)*(N2-1)+1);
		end
		if N2 == 1
			Dfull = zeros(Ne*(N1-1)+1,Ne);
		else
			Dfull = zeros(Ne*(N1-1)+1,(Ne)*(N2-1)+1);
		end
		Bfull = zeros(Ne*(N1-1)+1,2);

		for i = 1:Ne
			M1full[((i-1)*(N1-1)+1):((i)*(N1-1)+1),((i-1)*(N1-1)+1):((i)*(N1-1)+1)] = M1full[((i-1)*(N1-1)+1):((i)*(N1-1)+1),((i-1)*(N1-1)+1):((i)*(N1-1)+1)] + M1;
			L1full[((i-1)*(N1-1)+1):((i)*(N1-1)+1),((i-1)*(N1-1)+1):((i)*(N1-1)+1)] = L1full[((i-1)*(N1-1)+1):((i)*(N1-1)+1),((i-1)*(N1-1)+1):((i)*(N1-1)+1)] + L1(i);
			if N2 == 1
				M2full[i:i,i] = M2;
				L2full[i:i,i] = L2(i);
			else
				M2full[((i-1)*(N2-1)+1):((i)*(N2-1)+1),((i-1)*(N2-1)+1):((i)*(N2-1)+1)] = M2full[((i-1)*(N2-1)+1):((i)*(N2-1)+1),((i-1)*(N2-1)+1):((i)*(N2-1)+1)] + M2;
				L2full[((i-1)*(N2-1)+1):((i)*(N2-1)+1),((i-1)*(N2-1)+1):((i)*(N2-1)+1)] = L2full[((i-1)*(N2-1)+1):((i)*(N2-1)+1),((i-1)*(N2-1)+1):((i)*(N2-1)+1)] + L2(i);
			end
			if N2 == 1
				Dfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),i] = D;
			else
				Dfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),((i-1)*(N2-1)+1):((i)*(N2-1)+1)] = Dfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),((i-1)*(N2-1)+1):((i)*(N2-1)+1)]+D;
			end
			if i == 1
				Bfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),1] =Bfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),1]+B[:,1];
			end
			if i == Ne
				Bfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),2] =Bfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),2]+B[:,2];
			end
		end

		# interconnection matrix J:
		if N2 == 1
		phJ = [zeros(Ne*(N1-1)+1,Ne*(N1-1)+1) Dfull;
			-Dfull' zeros(Ne,Ne)];
		else
			phJ = [zeros(Ne*(N1-1)+1,Ne*(N1-1)+1) Dfull;
			-Dfull' zeros((Ne)*(N2-1)+1,(Ne)*(N2-1)+1)];
		end
		phD = zeros(2,2)
		# Hamiltonian matrix Q:
		Q = blkdiag((M1full\L1full)/(M1full), (M2full\L2full)/(M2full));

		# input matrix
		if N2 == 1
			phB = [Bfull; zeros(Ne,2)];
		else
			phB = [Bfull; zeros((Ne)*(N2-1)+1,2)];
		end

		p = Phs(phJ, phB, phD, Q)

		p
end

function weak_phs_FEM(Ne,N1, N2,a,b)
	# this function uses the partioned FEM method, finding a single PHS model
	# for the whole (1D) domain.
	# TO DO: - correct the input matrix, and allow for higher-order polynomials
	#	N1  (old = Norder+1)  variables x1, e1 and v1
	#	N2  (old = Norder)  #    variables x2, e2 and v2

		x1,w1,P = lglnodes(N1-1,(b-a)/Ne,0)    # discretization of x1 variables
		x1i,w1i,Pi = lglnodes(N1,0,(b-a)/Ne)  # these points are used to integrate the mass matrix using quadrature
		if N2 > 1
			x2,w2 = lglnodes(N2-1,(b-a)/Ne,0)
		else
			x2,w2 = [(a+b)/2], [(b-a)/Ne];
		end
		#M = massmatrix(ei,xi,we);
		#xquad, wquad = lgwt(Ne,a,b)\
		#M = massmatrix(ei,xi, xquad, wquad);
		M1 = massmatrix(x1,x1, x1i,w1i);
		if N2 == 1
			M2 = massmatrix(x2,x2, w2);
		else
			x2i,w2i,P2i = lglnodes(N2,0,(b-a)/Ne)  # these points are used to integrate the mass matrix using quadrature
			M2 = massmatrix(x2,x2, x2i,w2i);
		end

		D = transpose(dermatrix(x1,x2,1))*M2 ;

		p0 = map(i->leg_pol(0,x1,i), 1:length(x1))
		pL = map(i->leg_pol((b-a)/Ne,x1,i), 1:length(x1))
		B = [p0 pL]

		M1full = zeros(Ne*(N1-1)+1,Ne*(N1-1)+1);
		if N2 == 1
			M2full = zeros(Ne,Ne);
		else
			M2full = zeros((Ne)*(N2-1)+1,(Ne)*(N2-1)+1);
		end
		if N2 == 1
			Dfull = zeros(Ne*(N1-1)+1,Ne);
		else
			Dfull = zeros(Ne*(N1-1)+1,(Ne)*(N2-1)+1);
		end
		Bfull = zeros(Ne*(N1-1)+1,2);

		for i = 1:Ne
			M1full[((i-1)*(N1-1)+1):((i)*(N1-1)+1),((i-1)*(N1-1)+1):((i)*(N1-1)+1)] = M1full[((i-1)*(N1-1)+1):((i)*(N1-1)+1),((i-1)*(N1-1)+1):((i)*(N1-1)+1)] + M1;
			if N2 == 1
				M2full[i:i,i] = M2;
			else
				M2full[((i-1)*(N2-1)+1):((i)*(N2-1)+1),((i-1)*(N2-1)+1):((i)*(N2-1)+1)] = M2full[((i-1)*(N2-1)+1):((i)*(N2-1)+1),((i-1)*(N2-1)+1):((i)*(N2-1)+1)] + M2;
			end
			if N2 == 1
				Dfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),i] = D;
			else
				Dfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),((i-1)*(N2-1)+1):((i)*(N2-1)+1)] = Dfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),((i-1)*(N2-1)+1):((i)*(N2-1)+1)]+D;
			end
			if i == 1
				Bfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),1] =Bfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),1]+B[:,1];
			end
			if i == Ne
				Bfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),2] =Bfull[((i-1)*(N1-1)+1):((i)*(N1-1)+1),2]+B[:,2];
			end
		end

		# interconnection matrix J:
		if N2 == 1
		phJ = [zeros(Ne*(N1-1)+1,Ne*(N1-1)+1) Dfull;
			-transpose(Dfull) zeros(Ne,Ne)];
		else
			phJ = [zeros(Ne*(N1-1)+1,Ne*(N1-1)+1) Dfull;
			-transpose(Dfull) zeros((Ne)*(N2-1)+1,(Ne)*(N2-1)+1)];
		end
		phD = zeros(2,2)
		# Hamiltonian matrix Q:
		Q = blkdiag(inv(M1full), inv(M2full));

		# input matrix
		if N2 == 1
			phB = [Bfull; zeros(Ne,2)];
		else
			phB = [Bfull; zeros((Ne)*(N2-1)+1,2)];
		end

		p = Phs(phJ, phB, phD, Q)

		p
	end

function weak_phs_FEM(Ne,Norder,a,b)
	weak_phs_FEM(Ne,Norder+1,Norder+1,a,b)
end

#=
function weak_phs_FEM(Ne,Norder,a,b)
# this function uses the partioned FEM method, finding a single PHS model
# for the whole (1D) domain.
# TO DO: - correct the input matrix, and allow for higher-order polynomials
	N1 = Norder+1   # variables x1, e1 and v1
	N2 = Norder   	#    variables x2, e2 and v2

	x1,w1,P = lglnodes(N1-1,(b-a)/Ne,0)    # discretization of x1 variables
	x1i,w1i,Pi = lglnodes(N1,0,(b-a)/Ne)  # these points are used to integrate the mass matrix using quadrature
	if N2 > 1
		x2,w2 = lglnodes(N2-1,(b-a)/Ne,0)
	else
		x2,w2 = [(a+b)/2], [(b-a)/Ne];
	end
	#M = massmatrix(ei,xi,we);
	#xquad, wquad = lgwt(Ne,a,b)\
	#M = massmatrix(ei,xi, xquad, wquad);
	M1 = massmatrix(x1,x1, x1i,w1i);
	if Norder == 1
		M2 = massmatrix(x2,x2, w2);
	else
		x2i,w2i,P2i = lglnodes(N2,0,(b-a)/Ne)  # these points are used to integrate the mass matrix using quadrature
		M2 = massmatrix(x2,x2, x2i,w2i);
	end

	D = dermatrix(x1,x2,1)'*M2 ;

	p0 = map(i->leg_pol(0,x1,i), 1:length(x1))
	pL = map(i->leg_pol((b-a)/Ne,x1,i), 1:length(x1))
	B = [p0 pL]
	display(B)

	M1full = zeros(Ne*Norder+1,Ne*Norder+1);
	if Norder == 1
		M2full = zeros(Ne,Ne);
	else
		M2full = zeros((Ne)*(Norder-1)+1,(Ne)*(Norder-1)+1);
	end
	if Norder == 1
		Dfull = zeros(Ne*Norder+1,Ne);
	else
		Dfull = zeros(Ne*Norder+1,(Ne)*(Norder-1)+1);
	end
	Bfull = zeros(Ne*Norder+1,2);

	for i = 1:Ne
		M1full[((i-1)*Norder+1):((i)*Norder+1),((i-1)*Norder+1):((i)*Norder+1)] = M1full[((i-1)*Norder+1):((i)*Norder+1),((i-1)*Norder+1):((i)*Norder+1)] + M1;
		if Norder == 1
			M2full[i:i,i] = M2;
		else
			M2full[((i-1)*(Norder-1)+1):((i)*(Norder-1)+1),((i-1)*(Norder-1)+1):((i)*(Norder-1)+1)] = M2full[((i-1)*(Norder-1)+1):((i)*(Norder-1)+1),((i-1)*(Norder-1)+1):((i)*(Norder-1)+1)] + M2;
		end
		if Norder == 1
			Dfull[((i-1)*Norder+1):((i)*Norder+1),i] = D;
		else
			Dfull[((i-1)*Norder+1):((i)*Norder+1),((i-1)*(Norder-1)+1):((i)*(Norder-1)+1)] = Dfull[((i-1)*Norder+1):((i)*Norder+1),((i-1)*(Norder-1)+1):((i)*(Norder-1)+1)]+D;
		end
		if i == 1
			Bfull[((i-1)*Norder+1):((i)*Norder+1),1] =Bfull[((i-1)*Norder+1):((i)*Norder+1),1]+B[:,1];
		end
		if i == Ne
			Bfull[((i-1)*Norder+1):((i)*Norder+1),2] =Bfull[((i-1)*Norder+1):((i)*Norder+1),2]+B[:,2];
		end
	end

	# interconnection matrix J:
	if Norder == 1
	phJ = [zeros(Ne+1,Ne+1) Dfull;
		-Dfull' zeros(Ne,Ne)];
	else
		phJ = [zeros(Ne*Norder+1,Ne*Norder+1) Dfull;
		-Dfull' zeros((Ne)*(Norder-1)+1,(Ne)*(Norder-1)+1)];
	end
	phD = zeros(2,2)
	# Hamiltonian matrix Q:
	Q = blkdiag(inv(M1full), inv(M2full));

	# input matrix
	if Norder == 1
		phB = [Bfull; zeros(Ne,2)];
	else
		phB = [Bfull; zeros((Ne)*(Norder-1)+1,2)];
	end

	p = Phs(phJ, phB, phD, Q)

	p
end =#

function weak_phs1(N,a,b)
# this function uses the partioned FEM method, finding a PHS model for a
# part of the domain. This could be used with higher order polynomials
# (using a single element), or coupling several elements (using lower-order polynomials)
	N1 = N+1   # variables x1, e1 and v1
	N2 = N   	#    variables x2, e2 and v2

	x1,w1,P = lglnodes(N1-1,a,b)    # discretization of x1 variables
	if N2 > 1
		x2,w2 = lgwt(N2,a,b)
	else
		x2,w2 = (a+b)/2, (b-a);
	end
	#M = massmatrix(ei,xi,we);
	#xquad, wquad = lgwt(Ne,a,b)\
	#M = massmatrix(ei,xi, xquad, wquad);
	M1 = massmatrix(x1,x1, w1);
	M2 = massmatrix(x2,x2, w2);

	D = transpose(dermatrix(x1,x2,1))*M2 ;
	p0 = map(i->leg_pol(a,x1,i), 1:length(x1))
	pL = map(i->leg_pol(b,x1,i), 1:length(x1))

    phJ = [zeros(N1,N1) D;
	   -transpose(D) zeros(N2, N2)];
    phB = [p0 -pL; zeros(N2,2)];
    phD = zeros(2,2);

	Q = blkdiag(inv(M1),inv(M2));
	#pflow = Collocation_points(x1, w1)
	#peffort = Collocation_points(x2, w2)

	p = Phs(phJ, phB, phD, Q)
	#p.Q = Q
	#p.disc_data = Discrete_domain(pflow, peffort)
	#func(x :: Array) = (0.5*x'*p.Q*x)[1]
	#p.Hamiltonian = func
	#p.GradHam = ForwardDiff.gradient(p.Hamiltonian)
	p
end

function weak_phs2(N,a,b)
	N1 = N   # variables x1, e1 and v1
	N2 = N+1   	#    variables x2, e2 and v2

	x2,w2,P = lglnodes(N2-1,a,b)    # discretization of x1 variables
	if N1 > 1
		x1,w1 = lgwt(N1,a,b)
	else
		x1,w1 = (a+b)/2, (b-a);
	end
	#M = massmatrix(ei,xi,we);
	#xquad, wquad = lgwt(Ne,a,b)\
	#M = massmatrix(ei,xi, xquad, wquad);
	M1 = massmatrix(x1,x1, w1);
	M2 = massmatrix(x2,x2, w2);

	D = transpose(dermatrix(x2,x1,1))*M1 ;
	p0 = map(i->leg_pol(a,x2,i), 1:length(x2))
	pL = map(i->leg_pol(b,x2,i), 1:length(x2))

    phJ = [zeros(N1,N1) -transpose(D);
	   D zeros(N2, N2)];
    phB = [zeros(N1,2); p0 -pL];
    phD = zeros(2,2);

		Q = blkdiag(inv(M1),inv(M2));

	#pflow = Collocation_points(x1, w1)
	#peffort = Collocation_points(x2, w2)

	p = Phs(phJ, phB, phD, Q)
	#p.Q = Q
	#p.disc_data = Discrete_domain(pflow, peffort)
	#func(x :: Array) = (0.5*x'*p.Q*x)[1]
	#p.Hamiltonian = func
	#p.GradHam = ForwardDiff.gradient(p.Hamiltonian)
	p
end

#coupled_gyrator(ph1 :: Phs, ports1, ph2 :: Phs, ports2, couple_matrix)
function finelem(Nelem,Nint,a,b)
	dx = (b-a)/Nelem;
	p = weak_phs1(Nint,0,dx)
	for i = 2:Nelem
		if (i%2 == 1)
			pn = weak_phs1(Nint,0,dx)
		else
			pn = weak_phs2(Nint,0,dx)
		end
		p = PortHamiltonian.coupled_gyrator(p,[2],pn,[1],[-1])
	end
	p
end

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
module ProjSplx

export projsplx!, projsplx, projnorm1,
       projnorm1!, projnorm1t, proxinf!

"""
Projection onto the unit simplex {x | sum(x) = τ, x ≥ 0}.

   projsplx!(b, τ)

In-place variant of `projsplx`.
"""
function projsplx!{T}(b::Vector{T}, τ::T)

    n = length(b)
    bget = false

    idx = sortperm(b, rev=true)
    tsum = zero(T)

    @inbounds for i = 1:n-1
        tsum += b[idx[i]]
        tmax = (tsum - τ)/i
        if tmax ≥ b[idx[i+1]]
            bget = true
            break
        end
    end

    if !bget
        tmax = (tsum + b[idx[n]] - τ) / n
    end

    @inbounds for i = 1:n
        b[i] = max(b[i] - tmax, 0)
    end

end

"""
Projection onto the weighted simplex.

    projsplx!(b, c, τ)

This projects b (in-place) onto the unit simplex weighted by c.
"""
function projsplx!{T}(b::Vector{T}, c::Vector{T}, τ::T)

    n = length(b)
    bget = false

    @assert length(b) == length(c) "lengths must match"
    @assert minimum(c) > 0 "c is not positive."

    idx = sortperm(b./c, rev=true)
    tsum = csum = zero(T)

    @inbounds for i = 1:n-1
        j = idx[i]
        tsum += b[j]*c[j]
        csum += c[j]*c[j]
        tmax = (tsum - τ) / csum
        if tmax >= b[idx[i+1]] / c[idx[i+1]]
            bget = true
            break
        end
    end

    if !bget
        p = idx[n]
        tsum += b[p]*c[p]
        csum += c[p]*c[p]
        tmax = (tsum - τ) / csum
    end

    for i = 1:n
        @inbounds b[i] = max(b[i] - c[i]*tmax, 0)
    end

    return

end

"""
Projection onto the simplex.

  projsplx(b, τ) -> x

Variant of `projsplx`.
"""
function projsplx(b::Vector, τ)
    x = copy(b)
    projsplx!(x, τ)
    return x
end

"""
Projection onto the weighted simplex.

    projsplx(b, c, τ) -> x

Variant of `projsplx!`.
"""
function projsplx(b::Vector, c::Vector, τ)
    x = copy(b)
    projsplx!(x, c, τ)
    return x
end

# s = sign_abs!(x) returns
#     s[i] = true  if x[i] > 0
#     s[i] = false otherwise
# and x = abs(x).
function sign_abs!(x::Vector)
  n = length(x)
  s = Array{Bool}(n)
  @inbounds for i=1:n
    s[i] = x[i] ≥ 0
    x[i] = abs(x[i])
  end
  return s
end

# set_sign!(x, s) sets the sign of x based on s.
function set_sign!(x::Vector, s::Vector{Bool})
  n = length(x)
  @inbounds for i=1:n
      x[i] = s[i] ? x[i] : -x[i]
  end
end

"""
Projection onto the 1-norm ball

    {x | ||x||₁ ≤ τ }.

In-place variant of `projnorm1`.
"""
function projnorm1!(b::Vector, τ::Real)
    norm(b,1) > τ || return
    s = sign_abs!(b)
    projsplx!(b, τ)
    set_sign!(b, s)
end

"""
Projection onto the weighted 1-norm ball

    {x | ||diag(c)x||₁ ≤ τ }, c > 0.

In-place variant of `projnorm1`.
"""
function projnorm1!(b::Vector, c::Vector, τ::Real)
    norm(b,1) > τ || return
    s = sign_abs!(b)
    projsplx!(b, c, τ)
    set_sign!(b, s)
end

"""
Projection onto the weighted 1-norm ball.

Variant of `projnorm1!`.
"""
function projnorm1(b::Vector, c::Vector, τ::Real)
    x = copy(b)
    projnorm1!(x, c, τ)
    return x
end

"""
Projection onto the 1-norm ball.

Variant of `projnorm1!`.
"""
function projnorm1(b::Vector, τ::Real)
    x = copy(b)
    projnorm1!(x, τ)
    return x
end

"""
Proximal map of the scaled infinity norm.
    prox_inf(x,λ) = x - proj(x | λ𝔹₁)
     env_inf(x,λ) = (1/2λ)||x||² - (1/2λ)dist²(x | λ𝔹₁)
Modifies `x` in place; returns the envelope.
"""
function proxinf!(x::Vector, λ::Real)
  λ == 0 && return norm(x, Inf)
  nrmx2 = dot(x,x)
  xp = projnorm1(x, λ)
  BLAS.axpy!(-1., xp, x) # x <- x - xp
  return nrmx2/(2λ) - dot(x,x)/(2λ)
end

"""
Proximal map of the scaled infinity norm.
Return variant of `proxinf!`
"""
function proxinf(x::Vector, λ::Real)
    z = copy(x)
    proxinf!(z, λ)
    return z
end

end # module


text:
@testset "Invariant distribution tests" begin
   @testset "Embedding dim=$E" for E in 3:3
      @testset "Time series length: $ts_length" for ts_length in [10]
         @testset "Rep #$rep" for rep in 1:1
            println("E=",E, "\tlength(ts)=", ts_length, "\n")
            t_start = time_ns()

            # Create an invariant embedding
            embedding = gaussian_embedding(ts_length + E + 10)

            t = triang_from_embedding(embedding)

            # Split the triangulation such that all simplices have radii below some fraction
            # of the largest simplex radius in the triangulation.
            #print("Splitting triangulation...\n")
            #print("Initial number of simplices = ", size(t.simplex_inds, 1), "\n")

            #SimplexSplitting.refine_variable_k!(t, maximum(t.radii) - (maximum(t.radii) - mean(t.radii))/2)
            #print("Final number of simplices = ", size(t.simplex_inds, 1), "\n")

            #print("The number of simplices increased from...\n")

            # Create Markov matrix.
            #print("Testing sparse markov matrix computation\n")
            #t_sparse = time_ns()
            #M = mm_sparse(t)
            #print("mm_sparse took ", (time_ns() - t_sparse)/10^9, " seconds\n\n")
            print("Testing parallel markov matrix computation\n")
            t_parallel = time_ns()
            Mp = mm_p(t)
            print("mm_p took ", (time_ns() - t_parallel)/10^9, " seconds\n\n")

            t_markov = time_ns()

            # The column sums of the Markov matrix all must be 1. Otherwise, we haven't
            # accounted for the entire volume of the triangulation.
            @testset "Markov matrix" begin
               println(sum(Mp, 2))
               @test all(sum(Mp, 2) .≈ 1.0)
            end

            invmeasure, inds_nonzero_simplices = invariantdist(Array(Mp))
            t_measure = time_ns()
            #println("E=",E, "\tlength(ts)=", ts_length, "\trep=", rep, " | Markov: ", (t_markov - t_start)/10^9, " seconds for ", size(t.simplex_inds, 1), " simplices with ", nnz(M), " nonzero intersections, invmeasure: ", (t_measure - t_markov)/10^9, " seconds")

            # There must be at least one simplex with nonzero measure; otherwise, the
            # invariant distribution cannot sum to 1. Also, the invariant distribution must
            # sum to 1 within reasonable tolerance.
            @testset "Invariant measure" begin
               @test size(inds_nonzero_simplices, 1) >= 1
               @test sum(invmeasure) ≈ 1.0 atol = 1/10^6
            end
         end
      end
   end
end

max_stars_repo_name: kahaaga/InvariantDistribution

max_issues_repo_name: kahaaga/InvariantDistribution

max_forks_repo_name: kahaaga/InvariantDistribution


text:
function visualize(window::AbstractWindow, primitive)
    drawable, args = Drawable(window, primitive)
    push!(window[Renderlist], primitive => (drawable, args))
    return drawable
end


text:
module Callbacks

import PyCall: PyObject, pycall

import ..Keras
import ..Keras: PyDoc

const keras_callbacks = [
    "BaseLogger",
    "ProgbarLogger",
    "History",
    "ModelCheckpoint",
    "EarlyStopping",
    "RemoteMonitor",
    "LearningRateScheduler",
    "TensorBoard",
    "ReduceLROnPlateau",
    "CSVLogger",
    "LambdaCallback",
]

for c in keras_callbacks
    cb_name = Symbol(c)

    @eval begin
        struct $cb_name
            obj::PyObject

            @doc PyDoc(Keras._callbacks, Symbol($c)) function $cb_name(args...; kwargs...)
                new(Keras._callbacks[Symbol($c)](args...; kwargs...))
            end
        end

        PyObject(callback::$(cb_name)) = callback.obj
    end
end

end


text:
### INITIAL AND BOUNDARY CONDITIONS ###

#TODO: Retire DiscreteSpace and move to a DiscreteEquation class with corresponding DiscreteVariables which carry all the required information, define substitute methods for these.

abstract type AbstractBoundary end

abstract type AbstractTruncatingBoundary <: AbstractBoundary end

abstract type AbstractExtendingBoundary <: AbstractBoundary end

struct LowerBoundary <: AbstractTruncatingBoundary
    u
    x
    depvars
    indvars
    eq
    function LowerBoundary(u, t, x, eq, s, depvar_ops)
        depvars_lhs = get_depvars(eq.lhs, depvar_ops)
        depvars_rhs = get_depvars(eq.rhs, depvar_ops)
        depvars = collect(depvars_lhs ∪ depvars_rhs)
        #depvars =  filter(u -> !any(map(x-> x isa Number, arguments(u))), depvars)

        allx̄ = Set(filter(!isempty, map(u->filter(x-> t === nothing || !isequal(x, t.val), arguments(u)), depvars)))
        return new(u, x, depvar.(depvars, [s]), first(allx̄), eq)
    end
end

struct UpperBoundary <: AbstractTruncatingBoundary
    u
    x
    depvars
    indvars
    eq
    function UpperBoundary(u, t, x, eq, s, depvar_ops)
        depvars_lhs = get_depvars(eq.lhs, depvar_ops)
        depvars_rhs = get_depvars(eq.rhs, depvar_ops)
        depvars = collect(depvars_lhs ∪ depvars_rhs)

        allx̄ = Set(filter(!isempty, map(u->filter(x-> t === nothing || !isequal(x, t.val), arguments(u)), depvars)))
        return new(u, x, depvar.(depvars, [s]), first(allx̄), eq)
    end
end
struct PeriodicBoundary <: AbstractBoundary
    u
    x
end

getvars(b::AbstractBoundary) = (b.u, b.x)

@inline function isperiodic(bmps, u, x)
    if !isempty(bmps[operation(u)][x])
        # Being explicit hoping the compiler will optimize this away
        if first(bmps[operation(u)][x]) isa PeriodicBoundary
            return Val(true)
        else
            return Val(false)
        end
    else
        return Val(false)
    end
end

@inline function isperiodic(b)
    if b isa PeriodicBoundary
        return Val(true)
    else
        return Val(false)
    end
end

@inline function clip_interior!!(lower, upper, s, b::AbstractBoundary)
    # This x2i is correct
    dim = x2i(s, depvar(b.u, s), b.x)
    @assert dim !== nothing "Internal Error: Variable $(b.x) not found in $(depvar(b.u, s)), when parsing boundary condition $(b)"

    lower[dim] = lower[dim] + !isupper(b)
    upper[dim] = upper[dim] + isupper(b)
end

idx(b::LowerBoundary, s) = 1
idx(b::UpperBoundary, s) = length(s, b.x)

# indexes for Iedge depending on boundary type
isupper(::LowerBoundary) = false
isupper(::UpperBoundary) = true
isupper(::PeriodicBoundary) = false

@inline function edge(interiormap, s, u, j, islower)
    I = interiormap.I[interiormap.pde[depvar(u, s)]]
    # check needed on v1.6
    length(I) == 0 && return CartesianIndex{0}[]
    sd(i) = selectdim(I, j, i)
    I1 = unitindex(ndims(u, s), j)
    if islower
        edge = sd(1)
        # cast the edge of the interior to the edge of the boundary
        edge = edge .- [I1*(edge[1][j]-1)]
    else
        edge = sd(size(interiormap.I[interiormap.pde[depvar(u,s)]], j))
        edge = edge .+ [I1*(size(s.discvars[depvar(u, s)], j)-edge[1][j])]
    end
    return edge
end

edge(s, b, interiormap) = edge(interiormap, s, b.u, x2i(s, b.u, b.x), !isupper(b))

function _boundary_rules(s, orders, u, x, val)
        args = s.args[operation(u)]
        args = substitute.(args, (x=>val,))
        varrule = operation(u)(args...) => [operation(u)(args...), x]
        return vcat([(Differential(x)^d)(operation(u)(args...)) => [operation(u)(args...), x] for d in reverse(orders[x])], Differential(s.time)(operation(u)(args...)) => [operation(u)(args...), x], varrule)
end

function generate_boundary_matching_rules(s, orders)
    # TODO: Check for bc equations of multiple variables
    lowerboundary(x) = first(s.axies[x])
    upperboundary(x) = last(s.axies[x])

    # Rules to match boundary conditions on the lower boundaries
    lower = Dict([operation(u) => Dict([x => _boundary_rules(s, orders, u, x, lowerboundary(x)) for x in params(u, s)]) for u in s.ū])

    upper = Dict([operation(u) => Dict([x => _boundary_rules(s, orders, u, x, upperboundary(x)) for x in params(u, s)]) for u in s.ū])

    return (lower, upper)
end

"""
Creates a map of boundaries for each variable to be used later when discretizing the boundary condition equations, and
"""
function BoundaryHandler(bcs, s::DiscreteSpace, depvar_ops, tspan, derivweights::DifferentialDiscretizer)

    t=s.time

    if t === nothing
        initmaps = s.ū
    else
        initmaps = substitute.(s.ū,[t=>tspan[1]])
    end

    # Create some rules to match which bundary/variable a bc concerns
    # * Assume that the term of the condition is applied additively and has no multiplier/divisor/power etc.
    u0 = []
    bceqs = []
    ## BC matching rules, returns the variable and parameter the bc concerns

    lower_boundary_rules, upper_boundary_rules = generate_boundary_matching_rules(s, derivweights.orders)

    boundarymap = Dict([operation(u)=>Dict([x => [] for x in s.x̄]) for u in s.ū])

    # Generate initial conditions and bc equations
    for bc in bcs
        # * Assume in the form `u(...) ~ ...` for now
        depvarslhs = get_depvars(bc.lhs, depvar_ops)
        bcdepvar = first(depvarslhs)

        depvars = depvar.(collect(union(depvarslhs, get_depvars(bc.rhs, depvar_ops))), (s,))
        if any(u -> isequal(operation(u), operation(bcdepvar)), s.ū)
            if t !== nothing && ((operation(bc.lhs) isa Sym) | (operation(bc.lhs) isa Term)) && !any(x -> isequal(x, t.val), arguments(bc.lhs))
                # initial condition
                # * Assume that the initial condition is not in terms of the initial derivative i.e. equation is first order in time
                initindex = findfirst(isequal(bc.lhs), initmaps)
                if initindex !== nothing
                    #@show bcdepvar, bc, depvar(bcdepvar, s)
                    args = params(depvar(bcdepvar, s), s)
                    indexmap = Dict([args[i]=>i for i in 1:length(args)])
                    push!(u0,vec(s.discvars[s.ū[initindex]] .=> substitute.((bc.rhs,),valmaps(s, depvar(bcdepvar,s), depvars, indexmap))))
                end
            else
                # Split out additive terms
                terms = split_terms(bc)

                local u_, x_
                boundary = nothing
                # * Assume that the BC is defined on the edge of the domain
                # Check whether the bc is on the lower boundary, or periodic, we don't care which depvar/var
                for term in terms, r in reduce(vcat, reduce(vcat, collect.(values.(collect(values(lower_boundary_rules))))))
                    #Check if the rule changes the expression
                    if subsmatch(term, r)
                        # Get the matched variables from the rule
                        u_, x_ = r.second
                        # Mark the boundary
                        boundary = LowerBoundary(u_, s.time, x_, bc, s, depvar_ops)
                        # do it again for the upper end to check for periodic, but only check the current depvar and indvar
                        for term_ in setdiff(terms, [term]), r_ in upper_boundary_rules[operation(u_)][x_]
                            if subsmatch(term_, r_)
                                boundary = PeriodicBoundary(u_, x_)
                            end
                        end

                        break
                    end
                end
                # repeat for upper boundary
                if boundary === nothing
                    for term in terms, r in reduce(vcat, reduce(vcat, collect.(values.(collect(values(upper_boundary_rules))))))
                        if subsmatch(term, r)
                            u_, x_ = r.second
                            boundary = UpperBoundary(u_, s.time, x_, bc, s, depvar_ops)
                            break
                        end
                    end
                end
                @assert boundary !== nothing "Boundary condition $bc is not on a boundary of the domain, or is not a valid boundary condition"

                push!(boundarymap[operation(boundary.u)][boundary.x], boundary)
                #generate_bc_rules!(bceqs, Iedge, derivweights, s, bc, boundary)
            end
        end
    end
    return boundarymap, u0
end


text:
using Distributions
using ArgParse
# Import EnsembleKalmanProcesses modules
using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
using EnsembleKalmanProcesses.Observations
using EnsembleKalmanProcesses.ParameterDistributionStorage
include(joinpath(@__DIR__, "helper_funcs.jl"))

# Set parameter priors
param_names = ["cwtropics"]
n_param = length(param_names)
prior_dist = [Parameterized(Normal(35, 10))]
constraints = [[no_constraint()]]
priors = ParameterDistribution(prior_dist, constraints, param_names)

# Construct initial ensemble
N_ens = 10
initial_params = construct_initial_ensemble(priors, N_ens)
# Generate MiMA Parameters files
params_arr = [row[:] for row in eachrow(initial_params')]
versions = map(param -> generate_mima_params(param, param_names), params_arr)

# Store version identifiers for this ensemble in a common file
open("versions_1.txt", "w") do io
    for version in versions
        write(io, "mima_param_defs_$(version).txt\n")
    end
end


text:
using Documenter, BlockwiseInversion

makedocs(;
    modules=[BlockwiseInversion],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/djsegal/BlockwiseInversion.jl/blob/{commit}{path}#L{line}",
    sitename="BlockwiseInversion.jl",
    authors="djsegal",
    assets=[],
)

deploydocs(;
    repo="github.com/djsegal/BlockwiseInversion.jl",
)


text:
using Test
using TrajectoryOptimization
using RobotDynamics
using RobotDynamics: LieState
import RobotZoo.Quadrotor
using BenchmarkTools
using Rotations
using StaticArrays, LinearAlgebra, ForwardDiff
const TO = TrajectoryOptimization
const RD = RobotDynamics

import TrajectoryOptimization: DiagonalLieCost

## Test constructors
P = (3,6)
Qs = [rand(p) for p in P]
s = LieState(UnitQuaternion,3,6)
R = rand(4)
dcost = DiagonalLieCost(s, Qs, R)
@test dcost.w == [1]
@test dcost.Q == [Qs[1]; Qs[2]]
@test length(dcost.Q) == sum(P)

# Pass in Diagonal R
dcost = DiagonalLieCost(s, Qs, Diagonal(R))
@test dcost.R.diag == R

# Pass in Diagonal Q for vector parts
Q = Diagonal(vcat(Qs...))
q = rand(9)
dcost = DiagonalLieCost(s, Q, Diagonal(R), q=q, w=[5])
@test dcost.Q == [Qs[1]; Qs[2]]
@test dcost.q == q
@test dcost.w[1] === 5.0

# Pass in Diagonal Q for all states
Q = Diagonal(rand(13))
q = rand(13)
dcost = DiagonalLieCost(s, Q, Diagonal(R), q=q)
@test dcost.Q == [Q[1:3]; Q[8:13]]
@test dcost.q == [q[1:3]; q[8:13]]
@test dcost.w[1] ≈ sum(Q[4:7])

# Pass in weight, overriding default behavior of summing diagonal elements
dcost = DiagonalLieCost(s, Q, Diagonal(R), q=q, w=[6])
@test dcost.w[1] === 6.0

## Set up
model = Quadrotor{UnitQuaternion}()
x, u = rand(model)
s = LieState(model)
Q = rand(RobotDynamics.state_dim_vec(s))
R = rand(control_dim(model))
q = rand(RobotDynamics.state_dim_vec(s))
r = rand(control_dim(model))
c = rand()
w = rand(RobotDynamics.num_rotations(s))
qrefs = [rand(UnitQuaternion) for i = 1:RobotDynamics.num_rotations(s)]

## Call inner constructor
costfun = DiagonalLieCost(s, Q, R, q, r, c, w, qrefs)

# Test stage cost
TO.stage_cost(costfun, x)
p,quat,v,ω = RobotDynamics.parse_state(model, x)
Qr,Qv,Qω = Diagonal.((Q[1:3],Q[4:6],Q[7:9]))
Jv = 0.5*(p'Qr*p + v'Qv*v + ω'Qω*ω) + q[1:3]'p + q[4:6]'v + q[7:9]'ω
q0 = Rotations.params(quat)
qref = Rotations.params(qrefs[1])
dq = q0'qref
Jr = w[1]*min(1-dq,1+dq)
@test Jr + Jv ≈ TO.stage_cost(costfun, x)

Ju = 0.5*u'Diagonal(R)*u + r'u
@test TO.stage_cost(costfun, x, u) ≈ Jr + Jv + Ju

## Gradient
grad = ForwardDiff.gradient(x->TO.stage_cost(costfun, x), x)
gradq = -qref*w[1]*sign(dq)
@test grad ≈ [Qr*p + q[1:3]; gradq; Qv*v + q[4:6] ; Qω*ω + q[7:9]]
vinds = costfun.vinds
grad2 = zeros(length(grad))
grad2[vinds] .= Diagonal(Q)*x[vinds] + q

E = QuadraticCost{Float64}(13,4)
TO.gradient!(E, costfun, x, u)
@test E.q ≈ grad
@test E.r ≈ Diagonal(R)*u + r

## Hessian
hess = ForwardDiff.hessian(x->TO.stage_cost(costfun, x), x)
hess2 = zero(grad2)
hess2[vinds] .= Q
@test Diagonal(hess2) ≈ hess

E.Q .*= 0
@test TO.hessian!(E, costfun, x, u) == true
@test E.Q ≈ hess
@test E.R ≈ Diagonal(R)

## Try with an MRP
s1 = LieState(UnitQuaternion,P)
s2 = LieState(MRP,P)
dcost1 = DiagonalLieCost(s1, Qs, R)
dcost2 = DiagonalLieCost(s2, Qs, R)
model1 = Quadrotor{UnitQuaternion}()
model2 = Quadrotor{MRP}()
@test dcost2.Q == vcat(Qs...)
@test state_dim(dcost1) == 13
@test state_dim(dcost2) == 12
@test dcost2.vinds[end] == 12
@test dcost2.qinds[1] == [4,5,6]
u = @SVector rand(4)
x = rand(RBState)
x1 = RobotDynamics.build_state(model1, x)
x2 = RobotDynamics.build_state(model2, x)

@test TO.stage_cost(dcost1, x1, u) ≈ TO.stage_cost(dcost2, x2, u)

## Test conversion jacobians
q = rand(UnitQuaternion)

quat2mrp(q) = Rotations.params(MRP(UnitQuaternion(q)))
mrp2quat(g) = Rotations.params(UnitQuaternion(MRP(g)))
rp2quat(p) = Rotations.params(UnitQuaternion(RodriguesParam(p)))
quat2rp(q) = Rotations.params(RodriguesParam(UnitQuaternion(q)))

Rotations.jacobian(UnitQuaternion, MRP(q)) ≈
    ForwardDiff.jacobian(mrp2quat, Rotations.params(MRP(q)))
Rotations.jacobian(MRP, q) ≈
    ForwardDiff.jacobian(quat2mrp, Rotations.params(q))

Rotations.jacobian(UnitQuaternion, RodriguesParam(q)) ≈
    ForwardDiff.jacobian(rp2quat, Rotations.params(RodriguesParam(q)))

@test Rotations.jacobian(RodriguesParam, q) ≈
    ForwardDiff.jacobian(quat2rp, Rotations.params(q))

@btime Rotations.jacobian(UnitQuaternion, RodriguesParam($q))

g = Rotations.params(RodriguesParam(q))
rp2quat(g)
1/sqrt(1+g'g) * [1; g]

## Test solving a problem
model = Quadrotor()
N = 51
tf = 5.0
dt = tf / (N-1)

# Objective
x0,u0 = zeros(model)
xf = RD.build_state(model, [2,3,1], expm(SA[0,0,1]*deg2rad(135)), zeros(3), zeros(3))
Q = [fill(1,3), fill(0.1, 6)]
R = fill(1e-2,4)
costfun = TO.LieLQRCost(RD.LieState(model), Q, R, xf)
costfun_term = TO.LieLQRCost(RD.LieState(model), Q .* 100, R, xf)
obj = Objective(costfun, costfun_term, N)
obj2 = copy(obj)

prob = Problem(model, obj, xf, tf, x0=x0)
using Altro
solver = ALTROSolver(prob, show_summary=true)
solve!(solver)
states(solver)[end]
xf_sol = RBState(model, states(solver)[end])
x̄f = RBState(model, xf)
norm(xf_sol ⊖ x̄f)


text:
using Revise
using IMPES
using Plots

# using Plots, LinearAlgebra, IMPES
# using AnalyticalEOR

#* Define Relative Perms
kr_ww = SingleRelPerms(
        swr= 0.11,
        sor= 0.10,
        krw0= 0.2,
        kro0= 0.8,
        nw= 3.0,
        no= 2.0,
            )

kr_ow = SingleRelPerms(
        swr= 0.11,
        sor= 0.16,
        krw0= 0.4,
        kro0= 0.5,
        nw= 2.0,
        no= 3.0,
        )

#* Define reservoir
res = Reservoir(
    L = 300.,
    A = 200 * 50 * 0.3048^2,
    k = 100.0 * 9.869233e-16,
    phi = 0.2,
    )

#* Define fluids
fluid = Fluids(
            μw=1e-3,
            μo=5e-3,
)

#* Define grid
grid = Grid(
            res=res,
            Nx = 100,
            )

#* Define BC
bc = BoundaryConditions(
                        [:neumann, :dirichlet],
                        [426.5 * 0.3048 ^ 3 / 3600 / 24, 1000 * 6894.76]
                        )

tmax = 4000. * 3600 * 24 # days
dt = 1. * 3600 * 24

#* Initialize simulation and postprocessing
sim = Simulation(
                1000*ones(grid.Nx)* 6894.76,
                0.2*ones(grid.Nx),
                1*ones(grid.Nx),
                tmax,
                dt
)

post = FlowResults(sim, grid, tmax, dt)

while sim.t < tmax
    IMPES!(sim, res, kr_ww, fluid, grid, bc, post)
end

i = 1000
plot(grid.x, post.s[:,i], seriestype=:scatter)
plot(grid.x, post.p[:,i], seriestype=:scatter)

max_stars_repo_path: examples/Ex1 - BL1D.jl

max_issues_repo_path: examples/Ex1 - BL1D.jl

max_forks_repo_path: examples/Ex1 - BL1D.jl


text:
module smallPT

import Base: +, -, *, %
import Base.Threads: @threads
import ThreadPools: @qthreads

const inf = 1e20
const (DIFF, SPEC, REFR) = (1, 2, 3)

struct Vec
    x::Float64
    y::Float64
    z::Float64

    Vec(x = 0., y = 0. , z = 0.) = new(x, y, z)
end

const ZERO  = Vec(0., 0., 0.)
const XAXIS = Vec(1., 0., 0.)
const YAXIS = Vec(0., 1., 0.)

# vector operations

+(a::Vec, b::Vec) = return Vec(a.x + b.x, a.y + b.y, a.z + b.z)
-(a::Vec, b::Vec) = return Vec(a.x - b.x, a.y - b.y, a.z - b.z)
*(a::Vec, b::Float64) =  return Vec(a.x * b, a.y * b, a.z * b)
*(a::Float64, b::Vec) = b * a
*(a::Vec, b::Vec) = return Vec(a.x * b.x, a.y * b.y, a.z * b.z)
norm(a::Vec) = return a * (1 / sqrt(a.x^2 + a.y^2 + a.z^2))
dot(a::Vec, b::Vec) = return a.x * b.x + a.y * b.y + a.z * b.z                                      # dot product
%(a::Vec, b::Vec) = return Vec(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x) # cross product

struct Ray
    o::Vec
    d::Vec
end

struct Sphere
    rad::Float64
    p::Vec  # position
    e::Vec  # emission
    c::Vec  # color
    refl::Int   # reflection type: DIFFuse, SPECular , REFRactive
    # some precalculated values to speed up intersection computation
    maxRefl::Float64
    sqRad::Float64
    invC::Vec

    function Sphere(rad, p, e, c, refl)
        maxRefl = c.x > c.y && c.x > c.z ? c.x : c.y > c.z ? c.y : c.z    # max reflectivity
        sqRad = rad * rad
        new(rad, p, e, c, refl, maxRefl, sqRad, c*(1. / maxRefl) )
    end
end

function intersect(s::Sphere, r::Ray)
    op = s.p - r.o
    epsi = 1e-04
    b = dot(op, r.d)
    det = b * b - dot(op, op) + s.sqRad
    if det < 0.
        return inf
    else
        det = sqrt(det)
    end
    t1 = b - det
    t2 = b + det
    if  t1 > epsi
        return t1
    elseif t2 > epsi
        return t2
    end
    return inf
end

function cornell_box()

    return [Sphere(1e5, Vec(1e5+1, 40.8, 81.6), Vec(), Vec(.75, .25, .25), DIFF),      # Left
            Sphere(1e5, Vec(-1e5+99, 40.8, 81.6), Vec(), Vec(.25, .25, .75), DIFF),    # Right
            Sphere(1e5, Vec(50., 40.8, 1e5), Vec(), Vec(.75, .75, .75), DIFF),         # Back
            Sphere(1e5, Vec(50., 40.8, -1e5+170), Vec(), Vec(), DIFF),                 # Front
            Sphere(1e5, Vec(50., 1e5, 81.6), Vec(), Vec(.75, .75, .75), DIFF),         # Bottom
            Sphere(1e5, Vec(50., -1e5+81.6, 81.6), Vec(), Vec(.75, .75, .75), DIFF),   # Top
            Sphere(16.5, Vec(27., 16.5, 47.), Vec(), Vec(1., 1., 1.) *.999, SPEC),     # Mirror
            Sphere(16.5, Vec(73., 16.5, 78.), Vec(), Vec(1., 1., 1.)*.999, REFR),      # Glass
            Sphere(600, Vec(50., 681.6 - .27, 81.6), Vec(12.,12.,12.), Vec(), DIFF)    # Lite
    ]
end

function sky_scene()
    Cen = Vec(50.,40.8,-860.);
    return [
        Sphere(1600., Vec(1.,0.,2.)* 3000. , Vec(1.,.9,.8)*1.2e1*1.56*2., Vec(), DIFF), # sun
        Sphere(1560., Vec(1.,0.,2.)*3500.,Vec(1.,.5,.05)*4.8e1*1.56*2., Vec(),  DIFF), # horizon sun2
        #Sphere(10000,Cen+Vec(0,0,-200), Vec(0.0627, 0.188, 0.569)*6e-2*8., Vec(.7,.7,1)*.25,  DIFF), # sky
        Sphere(10000.,Cen+Vec(0.,0.,-200.), Vec(0.00063842, 0.02001478, 0.28923243)*6e-2*8., Vec(.7,.7,1.)*.25,  DIFF), # sky

        Sphere(100000., Vec(50., -100000., 0.),  Vec(),Vec(.3,.3,.3),DIFF), # grnd
        Sphere(110000., Vec(50., -110048.5, 0.),  Vec(.9,.5,.05)*4.,Vec(),DIFF),# horizon brightener
        Sphere(4e4, Vec(50., -4e4-30, -3000.),  Vec(),Vec(.2,.2,.2),DIFF),# mountains
        #Sphere(3.99e4, Vec(50., -3.99e4+20.045, -3000.),  Vec(),Vec(.7,.7,.7),DIFF),# mountains snow

        Sphere(26.5,Vec(22.,26.5,42.),   Vec(),Vec(1.,1.,1.)*.596, SPEC), # white Mirr
        Sphere(13.,Vec(75.,13.,82.),   Vec(),Vec(.96,.96,.96)*.96, REFR),# Glas
        Sphere(22.,Vec(87.,22.,24.),   Vec(),Vec(.6,.6,.6)*.696, REFR)    # Glas2
    ]
end

function vista()
    Cen = Vec(50.,-20.,-860.);
    return [
        Sphere(8000., Cen+Vec(0.,-8000.,-900.),Vec(1.,.4,.1)*5e-1, Vec(),  DIFF), # sun
        Sphere(1e4,  Cen+Vec(), Vec(0.631, 0.753, 1.00)*3e-1, Vec(1.,1.,1.)*.5,  DIFF), # sky

        Sphere(150.,  Cen+Vec(-350.,0., -100.),Vec(),  Vec(1.,1.,1.)*.3,  DIFF),    # mnt
        Sphere(200.,  Cen+Vec(-210.,0.,-100.), Vec(),  Vec(1.,1.,1.)*.3,  DIFF),    # mnt
        Sphere(145.,  Cen+Vec(-210.,85.,-100.),Vec(),  Vec(1.,1.,1.)*.8,  DIFF),    # snow
        Sphere(150.,  Cen+Vec(-50.,0.,-100.),  Vec(),  Vec(1.,1.,1.)*.3,  DIFF),    # mnt
        Sphere(150.,  Cen+Vec(100.,0.,-100.),  Vec(),  Vec(1.,1.,1.)*.3,  DIFF),    # mnt
        Sphere(125.,  Cen+Vec(250.,0.,-100.),  Vec(),  Vec(1.,1.,1.)*.3,  DIFF),    # mnt
        Sphere(150.,  Cen+Vec(375.,0.,-100.),  Vec(),  Vec(1.,1.,1.)*.3,  DIFF),    # mnt

        Sphere(2500., Cen+Vec(0.,-2400.,-500.),Vec(),  Vec(1.,1.,1.)*.1,  DIFF),    # mnt base

        Sphere(8000., Cen+Vec(0.,-8000.,200.), Vec(),  Vec(.2,.2,1),    REFR),      # water
        Sphere(8000., Cen+Vec(0.,-8000.,1100.),Vec(),  Vec(0.,.3,0.),     DIFF),    # grass
        Sphere(8.   , Cen+Vec(-75., -5., 850.),Vec(),  Vec(0.,.3,0.),     DIFF),    # bush
        Sphere(30.,   Cen+Vec(0.,   23., 825.),Vec(),  Vec(1.,1.,1.)*.996, REFR),   # ball

        Sphere(30.,  Cen+Vec(200.,280.,-400.),  Vec(),  Vec(1.,1.,1.)*.8,  DIFF),   # clouds
        Sphere(37.,  Cen+Vec(237.,280.,-400.),  Vec(),  Vec(1.,1.,1.)*.8,  DIFF),   # clouds
        Sphere(28.,  Cen+Vec(267.,280.,-400.),  Vec(),  Vec(1.,1.,1.)*.8,  DIFF),   # clouds

        Sphere(40.,  Cen+Vec(150.,280.,-1000.),  Vec(),  Vec(1.,1.,1.)*.8,  DIFF),  # clouds
        Sphere(37.,  Cen+Vec(187.,280.,-1000.),  Vec(),  Vec(1.,1.,1.)*.8,  DIFF),  # clouds

        Sphere(40.,  Cen+Vec(600.,280.,-1100.),  Vec(),  Vec(1.,1.,1.)*.8,  DIFF),  # clouds
        Sphere(37.,  Cen+Vec(637.,280.,-1100.),  Vec(),  Vec(1.,1.,1.)*.8,  DIFF),  # clouds

        Sphere(37.,  Cen+Vec(-800.,280.,-1400.),  Vec(),  Vec(1.,1.,1.)*.8,  DIFF), # clouds
        Sphere(37.,  Cen+Vec(0.,280.,-1600.),  Vec(),  Vec(1.,1.,1.)*.8,  DIFF),    # clouds
        Sphere(37.,  Cen+Vec(537.,280.,-1800.),  Vec(),  Vec(1.,1.,1.)*.8,  DIFF),  # clouds
    ]
end

const spheres = cornell_box()

# Converts floats to integers to be saved in PPM file
cl(x::Float64) = return max(min(x,one(x)),zero(x))
toInt(x::Float64) = return floor(Int, cl(x)^(1. / 2.2) * 255 + .5)

function intersectSpheres(r::Ray)
    t = inf
    id = nothing
    for sphere in spheres
        d = intersect(sphere, r)
        if d < t
            t = d
            id = sphere
        end
    end
    return id, t
end

function radiance(r::Ray, depth::Int)
    obj,t = intersectSpheres(r)
    if isnothing(obj)
        return ZERO
    end

    # Russian roulette: stop the recursion randomly based on the surface reflectivity

    depth+=1
    isMaxDepth = depth > 100
    useRR = depth > 5
    isRR = useRR && rand() < obj.maxRefl

    if isMaxDepth || (useRR && ~isRR)
        return obj.e
    end

    x = r.o + r.d * t                               # ray intersection position
    n = norm(x-obj.p)                               # sphere normal
    nl = dot(n, r.d) < 0. ? n : n * -1.             # properly oriented surface normal
    f = (useRR && isRR) ? obj.invC : obj.c    # object color (BRDF modulator)

    if obj.refl == DIFF
        r1 = 2 * pi * rand()
        r2 = rand()
        r2s = sqrt(r2)
        w = nl
        wo = w.x < -0.1 || w.x > 0.1 ? YAXIS : XAXIS
        u = norm(wo % w)
        v = w % u
        d = norm(u * cos(r1) * r2s + v * sin(r1) * r2s + w * sqrt(1-r2))
        return obj.e + f * radiance(Ray(x, d), depth)
    elseif obj.refl == SPEC
        return obj.e + f * radiance(Ray(x, r.d - n * (2. * dot(n ,r.d))), depth)
    else # otherwise we have a dielectric (glass) surface
        reflRay = Ray(x, r.d - n * (2. * dot(n ,r.d)))  # ideal dielectric reflection
        into = dot(n, nl) > 0.                          # Ray from outside going in?
        nc = 1.
        nt = 1.5
        nnt = into ? nc/nt : nt / nc
        ddn = dot(r.d, nl)
        cos2t = 1 - nnt * nnt * (1 - ddn * ddn)

        # if total internal reflecttion, reflect
        if cos2t < 0
            return obj.e + f * radiance(reflRay, depth)
        else
            # otherwise, choose reflection or refraction
            tdir = norm(r.d * nnt - n * ((into ? 1 : -1) * (ddn * nnt + sqrt(cos2t))))
            a = nt - nc
            b = nt + nc
            R0 = (a * a) / (b * b)
            c = 1. - (into ? -ddn : dot(tdir,n))
            Re = R0 + (1 - R0) * c * c * c * c * c
            Tr = 1 - Re
            P = .25 + .5 * Re
            RP = Re / P
            TP = Tr / (1. - P)

            result = ZERO
            if depth > 2
                if rand() < P
                    result = radiance(reflRay, depth) * RP
                else
                    result = radiance(Ray(x, tdir), depth) * TP
                end
            else
                result = radiance(reflRay, depth) * Re + radiance(Ray(x, tdir), depth) * Tr
            end

            return obj.e + f * result
        end
    end
end

function main()
    w = 1024
    h = 768
    samps = 64
    cam = Ray(Vec(50., 52., 295.6), norm(Vec(0., -0.042612, -1.)))  # camera position, direction
    fov = .5135                                                     # Field of view angle
    cx = Vec(w * fov/h, 0., 0.)                                     # Horizontal (x) camera direction increment
    cy = norm(cx % cam.d) * fov                                     # Vector up for the camera , y direction increment
    c = Array{Vec}(undef, w * h)                                    # The image
    for i = 1:length(c)
        c[i] = Vec(0., 0., 0.)
    end

    @qthreads for y in 1:h  # Loop over image rows
#=         if Threads.threadid() == 1
            percent = round(100 * y/(h-1), digits=1)
            print(stderr, "\rRendering ($(samps*4)) $(percent)%")   # Print progress
        end =#
        for x in 1:w                                                    # Loop columns
            # For each pixel do 2x2 subsamples and "samps" samples per subsample
            for sy in 1:2
                i = (h - y) * w + x                                     # Calculate array index for pixel(x,y)
                for sx in 1:2
                    r = Vec()
                    for s in 1:samps
                        # Tent filter
                        r1 = 2 * rand()
                        r2 = 2 * rand()
                        dx = r1 < 1. ? sqrt(r1) - 1. : 1. - sqrt(2. - r1)
                        dy = r2 < 1. ? sqrt(r2) - 1. : 1. - sqrt(2. - r2)
                        # Compute ray direction
                        d = cx * (((sx-1 + .5 + dx) / 2 + x) / w - .5) +
                            cy * (((sy-1 + .5 + dy) / 2 + y) / h - .5) + cam.d
                        r = r + radiance(Ray(cam.o + d * 140., norm(d)), 0) * (1. / samps)
                    end
                    @inbounds c[i] += Vec(cl(r.x), cl(r.y), cl(r.z)) * .25
                end
            end
        end
    end

    open("image_julia.ppm", "w") do f
        print(f, "P3\n$w $h\n255\n")
        for i in 1:length(c)
            print(f, toInt(c[i].x))
            print(f, " ")
            print(f, toInt(c[i].y))
            print(f, " ")
            print(f, toInt(c[i].z))
            print(f, "\n")
        end
    end
end

end