text:
"""
Solves the incompressible Navier-Stokes equations in three dimensions using a
pseudo-spectral approach involving the FFT. Simulates the 3D Taylor-Green vortex
from vortex stretching to fully turbulent to dissipation. The approach follows
the "SpectralDNS in Python" paper by Mortensen and Langtangen:
https://arxiv.org/pdf/1602.03638v1.pdf

Momentum:           ∂u/∂t + (u ⋅ ∇) u = − 1/ρ ∇p + ν ∇²u

Incompressibility:  ∇ ⋅ u = 0

u :  Velocity (3d vector)
p :  Pressure
ν :  Kinematic Viscosity
ρ :  Density (here=1)
t :  Time
∇ :  Nabla operator (defining nonlinear convection, gradient and divergence)
∇²:  Laplace Operator

------

Scenario:

The flow is simulated on a Ω = (0, 2π)³ domain (in 3D). The initial condition is
according to

    u(t=0, x, y, z) =   sin(x)cos(y)cos(z)

    v(t=0, x, y, z) = − cos(x)sin(y)cos(z)

    w(t=0, x, y, z) = 0

No forcing is applied.

------

Expected Outcome:

(Following page 18 of
https://doc.nektar.info/tutorials/latest/incns/taylor-green-vortex/incns-taylor-green-vortex.pdf)

1. Laminar Vortex sheets

2. Vortex stretching

3. Rearrangement

4. Breakdown

5. Fully turbulent

6. Dissipation & Decay

------

Solution Strategy:

u  : The velocity field (3D vector) in spatial domain
λ  : The velocity field (3D vector) in Fourier domain

ω  : The vorticity field (3D vector) in spatial domain
η  : The vorticity field (3D vector) in Fourier domain

m  : The product of vorticity and velocity (3D vector) in spatial domain
ϕ  : The product of vorticity and velocity (3D vector) in Fourier domain

k  : The wavenumber vector (3D vector)

ψ  : The pressure in Fourier domain (1D scalar)

b  : The rhs of the ODE system in Fourier domain (3D vector)

i  : The imaginary unit (1D scalar)

0. Initialize the velocity vectors according to the IC and transform them into
   Fourier Domain

1. Compute the curl in Fourier Domain by means of a cross
   product with the wavenumber vector and imaginary unit

    η = ℱ(ω) = ℱ(∇ × u) = i k × λ

2. Transform the vorticity back to spatial domain

    ω = ℱ⁻¹(η)

3. Compute the "convection" by means of a cross product

    m = u × ω

4. Transform the "convection" back into Fourier domain

    ϕ = ℱ(m)

5. Perform a dealising on the convection in order to suppress
   unresolved wave numbers

6. Compute the (pseudo) "pressure" in Fourier domain

    ψ = (k ⋅ ϕ) / ||k||₂

7. Compute the rhs to the ODE system

    b = m - ν ||k||₂² λ - (k ψ) / ||k||₂

8. Advance the velocity in Fourier Domain by means of an
   Explicit Euler time step

   λ ← λ + Δt b

9. Transform the newly obtained velocity back into spatial
   domain

   u = ℱ⁻¹(λ)

10. [Optional] visualize the vorticity magnitude in spatial
   domain interactively

11. Repeat from (1.)

In total, it takes us three (three-dimensional) Fourier Transforms
per time iteration:
- Transformation of the curl to spatial domain
- Transformation of the "convection" to Fourier Domain
- Transformation of the velocity to spatial domain

It is called pseudo-spectral because some operations are performed
in Fourier Domain and some in the spatial domain.

-------

Stability Consideration

Choose the time step size carefully.

-----

Additional Hints:

Check out the Python Repo:
https://github.com/spectralDNS/spectralDNS

In particular this file:
https://github.com/spectralDNS/spectralDNS/blob/master/spectralDNS3D_short.py

"""

using FFTW
using GLMakie
using Statistics
using ProgressMeter

N_POINTS_P_AXIS = 50
KINEMATIC_VISCOSITY = 1.0 / 1_600
TIME_STEP_LENGTH = 0.02
N_TIME_STEPS = 1_700
PLOT_EVERY = 3

function cross_product!(
    res_x::Array{T, 3},
    res_y::Array{T, 3},
    res_z::Array{T, 3},
    a_x::Array{T, 3},
    a_y::Array{T, 3},
    a_z::Array{T, 3},
    b_x::Array{T, 3},
    b_y::Array{T, 3},
    b_z::Array{T, 3},
) where T<:Union{Float64, Complex{Float64}}
    res_x .= (
        a_y .* b_z
        .-
        a_z .* b_y
    )

    res_y .= (
        a_z .* b_x
        .-
        a_x .* b_z
    )

    res_z .= (
        a_x .* b_y
        .-
        a_y .* b_x
    )
end

function main()
    x_range = range(0.0, 2*pi, N_POINTS_P_AXIS)
    y_range = range(0.0, 2*pi, N_POINTS_P_AXIS)
    z_range = range(0.0, 2*pi, N_POINTS_P_AXIS)

    coordinates_x = [x for x in x_range, y in y_range, z in z_range]
    coordinates_y = [y for x in x_range, y in y_range, z in z_range]
    coordinates_z = [z for x in x_range, y in y_range, z in z_range]

    wavenumbers_1d_x = fftfreq(N_POINTS_P_AXIS) .* N_POINTS_P_AXIS
    wavenumbers_1d_y = fftfreq(N_POINTS_P_AXIS) .* N_POINTS_P_AXIS
    wavenumbers_1d_z = fftfreq(N_POINTS_P_AXIS) .* N_POINTS_P_AXIS

    wavenumbers_x = [k_x for k_x in wavenumbers_1d_x, k_y in wavenumbers_1d_y, k_z in wavenumbers_1d_z]
    wavenumbers_y = [k_y for k_x in wavenumbers_1d_x, k_y in wavenumbers_1d_y, k_z in wavenumbers_1d_z]
    wavenumbers_z = [k_z for k_x in wavenumbers_1d_x, k_y in wavenumbers_1d_y, k_z in wavenumbers_1d_z]

    wavenumbers_norm = sqrt.(
        wavenumbers_x.^2
        +
        wavenumbers_y.^2
        +
        wavenumbers_z.^2
    )

    wavenumbers_norm[iszero.(wavenumbers_norm)] .= 1.0
    normalized_wavenumbers_x = wavenumbers_x ./ wavenumbers_norm
    normalized_wavenumbers_y = wavenumbers_y ./ wavenumbers_norm
    normalized_wavenumbers_z = wavenumbers_z ./ wavenumbers_norm

    # Taylor-Green vortex initial condition
    velocity_x =   sin.(coordinates_x) .* cos.(coordinates_y) .* cos.(coordinates_z)
    velocity_y = - cos.(coordinates_x) .* sin.(coordinates_y) .* cos.(coordinates_z)
    velocity_z =   zeros((N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))

    # Pre-Plan the FFT for additional performance
    fft_operator = plan_fft(velocity_x; flags=FFTW.MEASURE)

    velocity_x =   sin.(coordinates_x) .* cos.(coordinates_y) .* cos.(coordinates_z)

    # Pre-Allocatearrays
    velocity_x_fft = fft_operator * velocity_x
    velocity_y_fft = fft_operator * velocity_y
    velocity_z_fft = fft_operator * velocity_z

    curl_x = zeros(Float64, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))
    curl_y = zeros(Float64, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))
    curl_z = zeros(Float64, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))

    curl_x_squared = zeros(Float64, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))
    curl_y_squared = zeros(Float64, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))
    curl_z_squared = zeros(Float64, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))

    curl_magnitude = zeros(Float64, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))

    curl_x_fft = zeros(Complex{Float64}, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))
    curl_y_fft = zeros(Complex{Float64}, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))
    curl_z_fft = zeros(Complex{Float64}, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))

    convection_x = zeros(Float64, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))
    convection_y = zeros(Float64, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))
    convection_z = zeros(Float64, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))

    convection_x_fft = zeros(Complex{Float64}, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))
    convection_y_fft = zeros(Complex{Float64}, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))
    convection_z_fft = zeros(Complex{Float64}, (N_POINTS_P_AXIS, N_POINTS_P_AXIS, N_POINTS_P_AXIS))

    # De-Aliasing is necessary to stabilize the simulation
    k_max_dealias = 2.0/3.0 * (N_POINTS_P_AXIS//2 + 1)
    dealias = (
        (abs.(wavenumbers_x) .< k_max_dealias)
        .*
        (abs.(wavenumbers_y) .< k_max_dealias)
        .*
        (abs.(wavenumbers_z) .< k_max_dealias)
    )

    # Prepare the visualization
    makie_plot_figure_window = Figure()
    display(makie_plot_figure_window)

    makie_plot_axis = Axis3(
        makie_plot_figure_window[1, 1],
        perspectiveness = 0.3,
    )

    makie_plot_3d_contour = contour!(
        makie_plot_axis,
        x_range,
        y_range,
        z_range,
        curl_magnitude,
        alpha = 0.05,
        levels = range(0.2, 1.0, length=4),
    )
    println("Hit enter when plot is ready")
    readline()

    @showprogress "Simulating & Animating ..." for t in 1:N_TIME_STEPS
        # (1) Compute the Curl in Fourier Domain
        cross_product!(
            curl_x_fft,
            curl_y_fft,
            curl_z_fft,
            im .* wavenumbers_x,
            im .* wavenumbers_y,
            im .* wavenumbers_z,
            velocity_x_fft,
            velocity_y_fft,
            velocity_z_fft,
        )

        # (2) Transform curl to spatial domain
        curl_x .= real(fft_operator \ curl_x_fft)
        curl_y .= real(fft_operator \ curl_y_fft)
        curl_z .= real(fft_operator \ curl_z_fft)

        # (3) Compute "Convection" in spatial domain
        cross_product!(
            convection_x,
            convection_y,
            convection_z,
            velocity_x,
            velocity_y,
            velocity_z,
            curl_x,
            curl_y,
            curl_z,
        )

        # (4) Transform "Convection" to Fourier Domain
        convection_x_fft .= fft_operator * convection_x
        convection_y_fft .= fft_operator * convection_y
        convection_z_fft .= fft_operator * convection_z

        # (5) Dealiasing on the higher wavenumbers
        convection_x_fft .*= dealias
        convection_y_fft .*= dealias
        convection_z_fft .*= dealias

        # (6) Compute the Pseudo-Pressure by a Divergence in Fourier Domain
        pressure_fft = (
            normalized_wavenumbers_x .* convection_x_fft
            +
            normalized_wavenumbers_y .* convection_y_fft
            +
            normalized_wavenumbers_z .* convection_z_fft
        )

        # (7) Assemble the rhs to the ODE system in Fourier Domain
        rhs_x_fft = (
            convection_x_fft
            -
            KINEMATIC_VISCOSITY * wavenumbers_norm.^2 .* velocity_x_fft
            -
            normalized_wavenumbers_x .* pressure_fft
        )
        rhs_y_fft = (
            convection_y_fft
            -
            KINEMATIC_VISCOSITY * wavenumbers_norm.^2 .* velocity_y_fft
            -
            normalized_wavenumbers_y .* pressure_fft
        )
        rhs_z_fft = (
            convection_z_fft
            -
            KINEMATIC_VISCOSITY * wavenumbers_norm.^2 .* velocity_z_fft
            -
            normalized_wavenumbers_z .* pressure_fft
        )

        # (8) Euler Step Update
        velocity_x_fft .+= rhs_x_fft * TIME_STEP_LENGTH
        velocity_y_fft .+= rhs_y_fft * TIME_STEP_LENGTH
        velocity_z_fft .+= rhs_z_fft * TIME_STEP_LENGTH

        # (9) Transform the velocities back to spatial domain
        velocity_x .= real(fft_operator \ velocity_x_fft)
        velocity_y .= real(fft_operator \ velocity_y_fft)
        velocity_z .= real(fft_operator \ velocity_z_fft)

        # VISUALIZATION
        if t % PLOT_EVERY == 0
            curl_x_squared .= curl_x.^2
            curl_y_squared .= curl_y.^2
            curl_z_squared .= curl_z.^2
            curl_magnitude .= sqrt.(
                curl_x_squared
                +
                curl_y_squared
                +
                curl_z_squared
            )

            delete!(makie_plot_3d_contour.parent, makie_plot_3d_contour)

            makie_plot_3d_contour = contour!(
                makie_plot_axis,
                x_range,
                y_range,
                z_range,
                curl_magnitude,
                alpha = 0.01,
                levels = range(1.0, 10.0, length=5),
            )

            sleep(0.02)

        end

    end

end

main()


text:
import MPI

function main()
    MPI.init()

    comm = MPI.COMM_WORLD

    println("Hello world, I am $(MPI.rank(comm)) of $(MPI.size(comm))")

    MPI.barrier(comm)

    MPI.finalize()
end

main()

ext: jl


text:
ENV["PYTHON_JL_RUNTIME_PYTHON"] = Sys.which("python")
using Pkg
using PyCall
using OrderedCollections
using LinearAlgebra
using SpikeSynchrony
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
SNN.@load_units

#include("../opt_single_cell_utils.jl")

import DataStructures
using JLD

global ngt_spikes
global opt_vec
global extremum
global extremum_param
global ngt_spikes
global fitness

using Plots
unicodeplots()

###
global cell_type="ADEXP"
global vecp=false
###
using SpikeNetOpt
SNO = SpikeNetOpt

(vmgtv,vmgtt,ngt_spikes,ground_spikes) = SNO.get_data()
println("Ground Truth")
plot(vmgtt[:],vmgtv[:]) |> display

#lower,upper=get_izhi_ranges()

lower,upper=get_adexp_ranges()
ALLEN_DURATION = 2000 * ms
ALLEN_DELAY = 1000 * ms

function loss(E,ngt_spikes,ground_spikes)
    spikes = raster_synchp(E)
    spikes = [s*ms for s in spikes]
    maxt = findmax(sort!(unique(vcat(spikes,ground_spikes))))[1]
    if size(spikes)[1]>1
        t, S = SPIKE_distance_profile(spikes, ground_spikes;t0=0,tf = maxt)
        spkd = SpikeSynchrony.trapezoid_integral(t, S)/(t[end]-t[1]) # == SPIKE_distance(y1, y2)
    else
        spkd = 1.0
    end
    delta = abs(size(spikes)[1] - ngt_spikes)
    return spkd+delta

end

function zz_(param)
    if cell_type=="IZHI"
        pp = SNN.IZParameter(;a = param[1], b = param[2], c = param[3], d = param[4])
        E = SNN.IZ(;N = 1, param = pp)
    end
    if cell_type=="ADEXP"
        adparam = SNN.ADEXParameter(;a = param[1],
            b = param[2],
            cm = param[3],
            v0 = param[4],
            τ_m = param[5],
            τ_w = param[6],
            θ = param[7],
            delta_T = param[8],
            v_reset = param[9],
            spike_delta = param[10])

        E = SNN.AD(;N = 1, param=adparam)
    end

    SNN.monitor(E, [:v])
    SNN.monitor(E, [:fire])
    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms
    current = current_search(cell_type,param,ngt_spikes)
    E.I = [current*nA]#[param[5]*nA]

    SNN.sim!([E]; dt =1*ms, delay=ALLEN_DELAY,stimulus_duration=ALLEN_DURATION,simulation_duration = ALLEN_DURATION+ALLEN_DELAY+443ms)
    vecplot(E, :v) |> display
    error = loss(E,ngt_spikes,ground_spikes)
    error
end

function z(x::AbstractVector)
    param = x
    params = SNN.IZParameter(;a = param[1], b = param[2], c = param[3], d = param[4])
    E = SNN.IZ(;N = 1, param = params)
    SNN.monitor(E, [:v])
    SNN.monitor(E, [:fire])
    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms
    E.I = [param[5]*nA]
    SNN.sim!([E]; dt =1*ms, delay=ALLEN_DELAY,stimulus_duration=ALLEN_DURATION,simulation_duration = ALLEN_DURATION+ALLEN_DELAY+500ms)
    error = loss(E,ngt_spikes)
    @show(error)
    error

end

function checkmodel(param)
    if cell_type=="IZHI"
        pp = SNN.IZParameter(;a = param[1], b = param[2], c = param[3], d = param[4])
        E = SNN.IZ(;N = 1, param = pp)
    end
    if cell_type=="ADEXP"

        adparam = SNN.ADEXParameter(;a = param[1],
            b = param[2],
            cm = param[3],
            v_rest = param[4],
            tau_m = param[5],
            tau_w = param[6],
            v_thresh = param[7],
            delta_T = param[8],
            v_reset = param[9],
            spike_height = param[10])
    end
    E = SNN.AD(;N = 1, param=adparam)

    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms
    E.I = [current_search(param,ngt_spikes)*nA]

    SNN.monitor(E, [:v])
    SNN.sim!([E]; dt =1*ms, delay=ALLEN_DELAY,stimulus_duration=ALLEN_DURATION,simulation_duration = ALLEN_DURATION+ALLEN_DELAY+443ms)
    if vecp
        vec = SNN.vecplot(E, :v)
        vec |> display
        vec
    end

end

ɛ = 0.125
options = GA(
    populationSize = 20,
    ɛ = ɛ,
    selection = ranklinear(1.5),#ss,
    crossover = intermediate(1.0),#line(1.0),#xovr,
    mutation = uniform(1.0),#domainrange(fill(1.0,ts)),#ms
)
@time result = Evolutionary.optimize(zz_,lower,upper, initd, options,
    Evolutionary.Options(iterations=20, successive_f_tol=25, show_trace=true, store_trace=true)
)
fitness = minimum(result)
println("GA: ɛ=$ɛ) => F: $(minimum(result))")# C: $(Evolutionary.iterations(result))")
extremum_param = Evolutionary.minimizer(result)
opt_vec = checkmodel(extremum_param)
println("probably jumbled extremum param")
println(extremum_param)
plot(opt_vec)|> display
plot(vmgtt[:],vmgtv[:])|> display
#, label = "randData", ylabel = "Y axis",color = :red, legend = :topleft, grid = :off, xlabel = "Numbers Rand")
#p = twiny()
#plot!(p,vec, label = "log(x)", legend = :topright, box = :on, grid = :off, xlabel = "Log values") |> display
#return fitness
#end
#end
#break
#end


text:
#! format: off

using SIIPExamples

using PowerSystems
using PowerSimulations
const PSI = PowerSimulations
using PowerSystemCaseBuilder
using Dates

using HiGHS #solver

solver = optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.5)

sys_DA = build_system(PSITestSystems, "modified_RTS_GMLC_DA_sys")

sys_RT = build_system(PSITestSystems, "modified_RTS_GMLC_RT_sys")

template_uc = template_unit_commitment()
set_device_model!(template_uc, ThermalStandard, ThermalStandardUnitCommitment)

template_ed = template_economic_dispatch(
    network = NetworkModel(StandardPTDFModel, PTDF = PTDF(sys_DA), use_slacks = true),#NetworkModel(CopperPlatePowerModel, use_slacks = true),
)

models = SimulationModels(
    decision_models = [
        DecisionModel(template_uc, sys_DA, optimizer = solver, name = "UC"),
        DecisionModel(template_ed, sys_RT, optimizer = solver, name = "ED"),
    ],
)

feedforward = Dict(
    "ED" => [
        SemiContinuousFeedforward(
            component_type = ThermalStandard,
            source = OnVariable,
            affected_values = [ActivePowerVariable],
        ),
    ],
)

DA_RT_sequence = SimulationSequence(
    models = models,
    ini_cond_chronology = InterProblemChronology(),
    feedforwards = feedforward,
)

sim = Simulation(
    name = "rts-test",
    steps = 2,
    models = models,
    sequence = DA_RT_sequence,
    simulation_folder = mktempdir(".", cleanup = true),
)

build!(sim)

execute!(sim, enable_progress_bar = false)

results = SimulationResults(sim);
uc_results = get_decision_problem_results(results, "UC"); # UC stage result metadata
ed_results = get_decision_problem_results(results, "ED"); # ED stage result metadata

read_variables(uc_results)

read_parameters(uc_results)

list_variable_names(uc_results)

list_parameter_names(uc_results)

Dict([
    v => read_variable(uc_results, v) for v in [
        "ActivePowerVariable__RenewableDispatch",
        "ActivePowerVariable__HydroDispatch",
        "StopVariable__ThermalStandard",
    ]
])

read_parameter(
    ed_results,
    "ActivePowerTimeSeriesParameter__RenewableFix",
    initial_time = DateTime("2020-01-01T06:00:00"),
    count = 5,
)

read_realized_variables(
    uc_results,
    ["ActivePowerVariable__ThermalStandard", "ActivePowerVariable__RenewableDispatch"],
)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: NREL-SIIP/Examples

max_issues_repo_name: NREL-SIIP/Examples

max_forks_repo_name: NREL-SIIP/Examples


text:
module LaplaceBeltramiCubedSphereTestsMPI
  using PartitionedArrays
  using Test
  using FillArrays
  using Gridap
  using GridapPETSc
  using GridapGeosciences

  include("../ConvergenceAnalysisTools.jl")
  include("../LaplaceBeltramiCubedSphereTests.jl")

  function petsc_options()
    """
    -ksp_type cg -ksp_rtol 1.0e-06 -ksp_atol 0.0
    -ksp_monitor -pc_type asm -sub_ksp_type preonly
    -sub_pc_type lu
    """
  end

  function main(parts)
    GridapPETSc.with(args=split(petsc_options())) do
       num_refs=[2,3,4,5]
       hs=[2.0/2^n for n in num_refs]
       model_args_series=zip(Fill(parts,length(num_refs)),num_refs)
       t = PartitionedArrays.PTimer(parts,verbose=true)
       PartitionedArrays.tic!(t)
       hs1,k1errors,s1=convergence_study(solve_laplace_beltrami,
                                         hs,
                                         model_args_series,1,8,PETScLinearSolver())
       PartitionedArrays.toc!(t,"LaplaceBeltramiCubedSphereConvergenceStudy")
       display(t)
       @test round(s1,digits=1) ≈ 1.0

       # Do garbage collection of all PETSc objects
       # set up during convergence_study
       GridapPETSc.gridap_petsc_gc()

       PartitionedArrays.tic!(t)
       hs2,k2errors,s2=convergence_study(solve_laplace_beltrami,
                                         hs,
                                         model_args_series,2,12,PETScLinearSolver())
       PartitionedArrays.toc!(t,"LaplaceBeltramiCubedSphereConvergenceStudy")
       display(t)
       @test round(s2,digits=1) ≈ 2.0
    end
  end
  prun(main,mpi,4)
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
export SolitaryWaveWhithamGreenNaghdi, SolitaryWGN

"""
    SolitaryWaveWhithamGreenNaghdi(param; kwargs)

Compute the Whitham-Green-Naghdi solitary wave with prescribed velocity.

# Arguments
- `param :: NamedTuple`: parameters of the problem containing velocity `c` \
and dimensionless parameters `ϵ` and `μ`, \
and mesh size `L` and number of collocation points `N`;
## Keywords (optional)
- `guess :: Vector{Real}`: initial guess for the surface deformation (if not provided, the exact formula for SGN is used);
- `x₀ :: Real`: center of solitary wave (if guess is not provided);
- `SGN :: Bool`: if `true` computes the Serre-Green-Naghdi (instead of Whitham-Green-Naghdi) solitary wave (consider `SolitaryWaveSerreGreenNaghdi` instead);
- `method :: Int`: equation used (between `1` and `4`);
- `iterative :: Bool`: inverts Jacobian through GMRES if `true`, LU decomposition if `false` (default is `false`);
- `verbose :: Bool`: prints numerical errors at each step if `true` (default is `false`);
- `max_iter :: Int`: maximum number of iterations of the Newton algorithm (default is `20`);
- `tol :: Real`: relative tolerance measured in ℓ∞ norm (default is `1e-10`);
- `ktol :: Real`: tolerance of the Krasny filter (default is `0`, i.e. no filtering);
- `gtol :: Real`: relative tolerance of the GMRES algorithm (default is `1e-10`);
- `dealias :: Int`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `q :: Real`: Newton algorithm modified with `u_{n+1}=q*u_{n+1}+(1-q)*u_n` (default is `1`);
- `α :: Real`: adds `α` times spectral projection onto the Kernel to the Jacobian (default is `0`).

# Return values
`(η,u,v,mesh)` with
- `η :: Vector{Float64}`: surface deformation;
- `u :: Vector{Float64}`: layer-averaged velocity;
- `v :: Vector{Float64}`: derivative of the velocity potential at the surface;
- `mesh :: Mesh`: mesh collocation points.

"""
function SolitaryWaveWhithamGreenNaghdi(
                param :: NamedTuple;
                guess = zeros(0) :: Vector{Float64},
                x₀ = 0 :: Real,
                SGN = false :: Bool,
                method = 2 :: Int,
                iterative = false :: Bool,
                verbose = false :: Bool,
                max_iter = 20 :: Int,
                tol = 1e-10 :: Real,
                ktol = 0 :: Real,
                gtol = 1e-10 :: Real,
                dealias = 0 :: Int,
                q=1 :: Real,
                α=0 :: Real
                        )

        c = param.c
        ϵ = param.ϵ
        μ = param.μ
        if abs(c)<1
                @error("The velocity must be greater than 1 (in absolute value).")
        end

        mesh = Mesh(param)

        if guess == zeros(0)
                # Using the exact formula for the SGN solitary wave as initial guess (if not provided)
                guess = (c^2-1)/ϵ*sech.(sqrt(3*(c^2-1)/(c^2)/μ)/2*(mesh.x.-x₀)).^2
        end

        k = mesh.k

        Dx       =  1im * k
        if SGN == true
                F₀ = √μ*Dx #./ (1 .+ μ/3 * k.^2 ).^(1/4)
        else
                F₁ 	= tanh.(sqrt(μ)*abs.(k))./(sqrt(μ)*abs.(k))
                F₁[1] 	= 1
                F₀       = 1im * sqrt.(3*(1 ./F₁ .- 1)).*sign.(k)
        end

        if dealias == 0
                Π = ones(size(k)) # no dealiasing (Π=Id)
        else
                K = (mesh.kmax-mesh.kmin)/(2+dealias)
                Π = abs.(k) .<= K # Dealiasing low-pass filter
        end

        krasny(k) = (abs.(k).> ktol ).*k
        krasny!(k) = k[abs.(k).< ktol ].=0

        function filter( v :: Vector{Float64} )
                real.(ifft(krasny(Π.*fft(v))))
        end
        function filter( v :: Vector{Complex{Float64}} )
                ifft(krasny(Π.*fft(v)))
        end
        function F( v , hv, Fv, F2v )
                if method == 1
                        return -1/3 ./ hv.^2 .* F2v .+
                                ϵ/2 .* (hv.*Fv).^2 .+
                                v./hv .- 1/(c^2) .*v .- ϵ/2 .* (v./hv).^2
                elseif method == 2
                        return -c/3 ./ hv.^2 .* F2v .+
                                ϵ/2 .* Fv.^2 .+
                                c*v .- 1/c .*v.*hv .- ϵ/2 .* v.^2
                elseif method == 3
                        return -1/3 ./ hv.^2 .* F2v .+
                                ϵ/2 .* Fv.^2 .+
                                v .- 1/(c^2) .*v.*hv .- ϵ/2 .* v.^2
                elseif method == 4
                        return (-1/3 .* F2v .+
                                ϵ/2 .*  hv.^4 .* Fv.^2 .+
                                hv.^2 .* (v .- ϵ/2 .* v.^2) .-
                                1/(c^2) .* hv.^3 .* v)/c^4
                end
        end

        function Fabs( v , hv, Fv, F2v )
                if method == 1
                        return abs.(1/3 ./ hv.^2 .* F2v) .+
                                abs.(ϵ/2 .* (hv.*Fv).^2 ) .+
                                abs.(v./hv) .- 1/(c^2) .*abs.(v) .+ ϵ/2 .* abs.((v./hv).^2)
                elseif method == 2
                        return abs.(c/3 ./ (hv.^2).* F2v) .+
                                abs.(ϵ/2 .* Fv.^2) .+
                                abs.(c*v) .+ abs.(1/c .*v.*hv) .+ abs.(ϵ/2 .* v.^2)
                elseif method == 3
                        return abs.(1/3 ./ (hv.^2).* F2v) .+
                                abs.(ϵ/2 .* Fv.^2) .+
                                abs.(v) .+ abs.(1/(c^2) .*v.*hv) .+ abs.(ϵ/2 .* v.^2)
                elseif method == 4
                        return (abs.(-1/3 .* F2v) .+
                                abs.(ϵ/2 .*  hv.^4 .* Fv.^2) .+
                                abs.( hv.^2 .*v) .+ abs.( ϵ/2 .* hv.^2 .* v.^2) .+
                                1/(c^2) .* abs.(hv.^3 .* v))/c^4
                end
        end
        if iterative == false
                x = mesh.x
                x₀ = mesh.x[1]
                FFT = exp.(-1im*k*(x.-x₀)');
                IFFT = exp.(1im*k*(x.-x₀)')/length(x);
                #Id = Diagonal(ones(size(x)));
                M₀ = IFFT * Diagonal( F₀ )* FFT
                M(v) = Diagonal( v )
                function JacF( v , hv, Fv , F2v, Dv, dxv )
                        if method == 1
                                return Symmetric(real.(
                                        -1/3 *M(1 ./ hv.^2)* M₀ * M(hv.^3)* ( M₀ * M( 1 ./ hv.^2 ) .+ 3* ϵ * M( Fv ./hv ) )
                                        .+ ϵ * M( hv.^2 .* Fv ) * M₀ * M( 1 ./ hv.^2 )
                                        .+ M( Dv ) .+ α*dxv*dxv' ))
                        elseif method == 2
                                return real.(
                                        -1/3 *M(1 ./ hv.^2 )* M₀ * M(hv.^3)* (c*M₀ .+ 3*ϵ * M( Fv ))
                                        .+ ϵ * M( hv .* Fv ) * M₀
                                        #.+ M(1 ./hv) .+ 2*ϵ/3 * M(F2v ./ hv) .+ ϵ^2 * M(hv .* Fv.^2) .- M((hv/c).^2)
                                        .+ M( Dv ) .+ α*M( hv.^2 )*dxv*dxv' *M(1 ./ hv.^2 ) )
                        elseif method == 3
                                return real.(
                                        -1/3 *M(1 ./ hv.^2 )* M₀ * M(hv.^3)* (M₀ .+ 3*ϵ * M( Fv ))
                                        .+ ϵ * M( hv .* Fv ) * M₀
                                        #.+ M(1 ./hv) .+ 2*ϵ/3 * M(F2v ./ hv) .+ ϵ^2 * M(hv .* Fv.^2) .- M((hv/c).^2)
                                        .+ M( Dv ) .+ α*M( hv.^2 )*dxv*dxv' *M(1 ./ hv.^2 ) )
                        elseif method == 4
                                return Symmetric(real.(
                                ( M₀ * M(hv.^3 ) * ( -1/6*M₀ .- ϵ * M( hv .* Fv ) ) .+
                                 ( -1/6*M₀ .+ ϵ * M( hv .* Fv ) ) * M( hv.^3  ) * M₀ .+
                                 M(hv.^2)*M( Dv ))/c^4 .+ α*dxv*dxv'
                                ))
                        end
                end
        else
                Four(v) = (ifft(F₀.*fft(v)))
                function JacFfast( v , hv, Fv , F2v, Dv, dxv )
                        function dF(φ)
                                if method == 1
                                 return fft( -1/3 ./ hv.^2 .* Four( hv.^3 .*
                                                ( Four( ifft(φ) ./ hv.^2 ) .+ 3* ϵ * ifft(φ) .* Fv ./hv )  )
                                        .+ ϵ .* hv.^2 .* Fv .* Four( ifft(φ) ./ hv.^2 )
                                         .+ ifft(φ) .* Dv .+ α*dot(dxv,ifft(φ))*dxv )
                                 elseif method == 2
                                  return fft( -1/3 ./ hv.^2 .* Four( hv.^3 .* ( c*Four( ifft(φ) ) .+ 3*ϵ*Fv.*ifft(φ) ))
                                          .+ ϵ .* hv .* Fv .* Four( ifft(φ) )
                                          .+ ifft(φ) .* Dv .+ α*dot(dxv, ifft(φ) ./ hv.^2)*hv.^2 .*dxv )
                                  elseif method == 3
                                   return fft( -1/3 ./ hv.^2 .* Four( hv.^3 .* ( Four( ifft(φ) ) .+ 3*ϵ*Fv.*ifft(φ) ))
                                           .+ ϵ .* hv .* Fv .* Four( ifft(φ) )
                                           .+ ifft(φ) .* Dv .+ α*dot(dxv, ifft(φ) ./ hv.^2)*hv.^2 .*dxv )
                                elseif method == 4
                                 return fft(( -1/3  .* Four( hv.^3 .* ( Four( ifft(φ) ) .+ 3*ϵ*Fv.*hv.*ifft(φ) ))
                                        .+ ϵ .* hv.^4 .* Fv .* Four( ifft(φ) )
                                        .+ ifft(φ) .* hv.^2 .* Dv)/c^4 .+ α*dot(dxv,ifft(φ))*dxv)
                                end
                         end
                if method == 1 || method == 4
                        return LinearMap(dF, length(v); issymmetric=true, ismutating=false)
                else
                        return LinearMap(dF, length(v); issymmetric=false, ismutating=false)
                end
                end

        end
#------ Iteration scheme

        # initial guess for iteration
        if method == 1
                u = filter(guess)
        elseif method == 2
                u = c*filter(guess./(1 .+ ϵ*guess))
        elseif method == 3 || method == 4
                u = filter(guess./(1 .+ ϵ*guess))
        end

        du,fu,dxu,hu,Fu,F2u,Du = (similar(u),).*ones(7)

        for i in range(0, stop=max_iter)
                dxu .= real.(ifft(Dx.*fft(u)))
                dxu ./= norm(dxu,2)
                if method == 1
                        hu .= 1 .+ϵ*u
                        Fu .= real.(ifft(F₀.*fft(u./hu)))
                        F2u .= real.(ifft(F₀.*fft( hu.^3 .* Fu ) ))
                        Du .= 1 ./ hu.^3 .+ 2*ϵ/3 ./ hu.^3 .* F2u .+ ϵ^2 * hu.* Fu.^2 .- 1/c^2
                elseif method == 2
                        hu .= c ./(c .- ϵ*u)
                        Fu .= hu.* real.(ifft(F₀.*fft(u)))
                        F2u .= real.(ifft(F₀.*fft(hu.^2 .* Fu ) ))
                        Du .= c ./hu .+ 2*ϵ/3 * F2u ./ hu .+ ϵ^2/c * hu .* Fu.^2 .- (hu.^2)/c
                elseif method == 3
                        hu .= 1 ./(1 .- ϵ*u)
                        Fu .= hu.* real.(ifft(F₀.*fft(u)))
                        F2u .= real.(ifft(F₀.*fft(hu.^2 .* Fu ) ))
                        Du .= 1 ./hu .+ 2*ϵ/3 * F2u ./ hu .+ ϵ^2 * hu .* Fu.^2 .- (hu/c).^2
                elseif method == 4
                        hu .= 1 ./(1 .- ϵ*u)
                        Fu .= real.(ifft(F₀.*fft(u)))
                        F2u .= real.(ifft(F₀.*fft(hu.^3 .* Fu ) ))
                        Du .=  1 ./ hu .+ 2*ϵ^2 * hu.^3 .* Fu.^2  .+ 2*ϵ.* hu .*(u.- ϵ/2* u.^2).-(hu .+ 3*ϵ* u.* hu.^2 )./(c^2)
                end
                fu .= F(u,hu,Fu,F2u)

                relerr = norm(fu,Inf)/norm(Fabs(u,hu,Fu,F2u),Inf)
                abserr = norm(fu,Inf)
    		if relerr < tol
    			@info string("Converged : relative error ",relerr," in ",i," steps\n")
    			break
    		elseif verbose == true
                        print(string("absolute error at step ",i,": ",abserr,"\n"))
                        print(string("relative error at step ",i,": ",relerr,"\n"))

    		end
                if i == max_iter
                        @warn string("The algorithm did not converge after ",i," steps: final relative error is ",relerr,"\n")
                        break
                end

                if iterative == false
                        du .=  JacF(u,hu,Fu,F2u,Du,dxu)  \ fu
                else
                        x = mesh.x
                        x₀ = mesh.x[1]
                        FFT = exp.(-1im*k*(x.-x₀)');
                        IFFT = exp.(1im*k*(x.-x₀)')/length(x);
                        F₁ 	= tanh.(sqrt(μ)*abs.(k))./(sqrt(μ)*abs.(k))
                        F₁[1] 	= 1
                        F₁ = 1 ./ (1 .+ μ*k.^2   )
                        #F₁ = ones(length(k))
                        Precond = Diagonal( 1 ./ F₁  )
                        du .=  real.(ifft(gmres( JacFfast(Complex.(u),hu,Fu,F2u,Du,dxu) , fft(fu) ; Pl = Precond, reltol = gtol, verbose=verbose )))
                end
    		u .-= q*filter(du)
        end

        if method == 1
                η = u
        elseif method ==2
                η = u./(c .- ϵ*u)
        elseif method ==3 || method == 4
                η = u./(1 .- ϵ*u)
        end

        h = 1 .+ ϵ*η
        u = c*η./h
        DxF(v) = real.(ifft(F₀ .* fft(v)))
	v = u - 1/3 ./h .* (DxF(h.^3 .*DxF(u)))

        return (η,u,v,mesh)

end

"""
    SolitaryWGN(param; kwargs)

Build the initial data associated with `SolitaryWaveWhithamGreenNaghdi(param; kwargs)`, of type `InitialData`,
to be used in initial-value problems `Problem(model, initial::InitialData, param)`.

---
    SolitaryWGN(c; ϵ=1,μ=1,N=2^12,kwargs)
Build the initial data with velocity `c`, dimensionless parameters `ϵ` and `μ`, and number of collocation points `N`, \
and `kwargs` the other (optional) keyword arguments as above.
"""
struct SolitaryWGN <: InitialData

	η
	v
	label :: String
	info  :: String

	function SolitaryWGN(param::NamedTuple;
                                guess = zeros(0) :: Vector{Float64},
                                x₀ = 0 :: Real,
                                SGN = false :: Bool,
                                method = 2 :: Int,
                                iterative = false :: Bool,
                                verbose = false :: Bool,
                                max_iter = 20 :: Int,
                                tol = 1e-10 :: Real,
                                ktol = 0 :: Real,
                                gtol = 1e-10 :: Real,
                                dealias = 0 :: Int,
                                q=1 :: Real,
                                α=0 :: Real)
		(η,u,v,mesh)=SolitaryWaveWhithamGreenNaghdi(param;
                                guess=guess,x₀=x₀,SGN=SGN,method=method,
                                iterative=iterative,verbose=verbose,max_iter=max_iter,
                                tol=tol,ktol=ktol,gtol=gtol,dealias=dealias,q=q,α=α)
		init = Init(mesh,η,v)
                if SGN == true
        	    model = "Serre-Green-Naghdi"
                else
                    model = "Whitham-Green-Naghdi"
                end
                label = "$model solitary wave"
		info = "Solitary travelling wave for the $model model.\n\
		├─velocity c = $(param.c)\n\
		└─maximum h₀ = $(maximum(η)) (from rest state)."

		new( init.η,init.v,label,info  )
	end

	function SolitaryWGN(
                        c::Real; ϵ=1::Real,μ=1::Real,N=2^12::Int,
                        guess = zeros(0) :: Vector{Float64},
                        x₀ = 0 :: Real,
                        SGN = false :: Bool,
                        method = 2 :: Int,
                        iterative = false :: Bool,
                        verbose = false :: Bool,
                        max_iter = 20 :: Int,
                        tol = 1e-10 :: Real,
                        ktol = 0 :: Real,
                        gtol = 1e-10 :: Real,
                        dealias = 0 :: Int,
                        q=1 :: Real,
                        α=0 :: Real
                                )
		L=200/sqrt(3*(c^2-1)/(c^2)/μ)/2
		xmin,xmax = x₀-L,x₀+L;
		param=(ϵ=ϵ,μ=μ,c=c,L=L,xmin=xmin,xmax=xmax,N=N)
		sol=SolitaryWGN(param;
                        guess=guess,x₀=x₀,SGN=SGN,method=method,
                        iterative=iterative,verbose=verbose,max_iter=max_iter,
                        tol=tol,ktol=ktol,gtol=gtol,dealias=dealias,q=q,α=α)
		new( sol.η,sol.v,sol.label,sol.info  )
	end

end


text:
using AngularMomentumAlgebra.Dipoles

@testset "Wigner–Eckart" begin
    @testset "Dipoles" begin
        x,y,z = 𝐫̂

        os = os"k[s-h]"

        @testset "z" begin
            for o in os
                ℓ = o.ℓ
                mℓs = -ℓ:ℓ
                for mℓ in mℓs
                    @test matrix_element((ℓ,mℓ),z,(ℓ+1,mℓ)) ≈ √((ℓ+mℓ+1)*(ℓ-mℓ+1)/((2ℓ+3)*(2ℓ+1)))
                end
            end
        end

        @testset "x" begin
            for o in os
                ℓ = o.ℓ
                mℓs = -ℓ:ℓ
                for mℓ in mℓs
                    a = SpinOrbital(o, mℓ, half(1))
                    b = SpinOrbital(Orbital(o.n, ℓ+1), mℓ+1, half(1))
                    # TODO Figure out where the minus sign comes from
                    @test dot(a,x,b) ≈ -0.5*√((ℓ+mℓ+2)*(ℓ+mℓ+1)/((2ℓ+3)*(2ℓ+1)))
                end
            end
        end
    end
end


text:
using Test
using DataTypesBasic

@test isempty(detect_ambiguities(DataTypesBasic))

@testset "Identity" begin
  include("Identity.jl")
end
@testset "Const" begin
  include("Const.jl")
end
@testset "Option" begin
  include("Option.jl")
end
@testset "Either" begin
  include("Either.jl")
end
@testset "Try" begin
  include("Try.jl")
end
@testset "ContextManager" begin
  include("ContextManager.jl")
end

@testset "convert" begin
  include("convert.jl")
end
@testset "promote_type" begin
  include("promote_type.jl")
end


text:
using NumericalDataManipulation.Data.Storage

@enum MergeFunctionStrategyName NAIVE_JOIN SIGMOID_JOIN

struct GridMergeParameters
    tail_knots::Int
    ϵʳᵉˡ::Float64
    meanΔ::Function
end

abstract type Task end
struct MergeTwoTablesTask <: Task
    file_a::FileResource
    file_b::FileResource
    master::Symbol
    columns::Vector{Int}
    master_interval::Tuple{Float64, Float64, Symbol}
    slave_intervals::Vector{Tuple{Float64, Float64, Symbol}}
    merge_function_strategy::MergeFunctionStrategyName
    grid_merge_parameters::GridMergeParameters
    skip_title_a::Bool
    skip_title_b::Bool
    function MergeTwoTablesTask(
        path_a::AbstractString, path_b::AbstractString, master::Symbol,
        columns::Vector{Int},
        master_interval::Tuple{Float64, Float64, Symbol},
        slave_intervals::Vector{Tuple{Float64, Float64, Symbol}},
        merge_function_strategy::MergeFunctionStrategyName,
        grid_merge_parameters::GridMergeParameters,
        skip_title_a::Bool, skip_title_b::Bool)
        return new(FileResource(path_a), FileResource(path_b), master,
            columns, master_interval, slave_intervals,
            merge_function_strategy, grid_merge_parameters,
            skip_title_a, skip_title_b)
    end
end


text:
using Test, JuLIP, JuLIPMaterials, LinearAlgebra
using JuLIPMaterials: ForceConstantMatrix1

(@isdefined test_pots) || (test_pots = JuLIP.Deps.fetch_test_pots())

eam_Fe = EAM(test_pots * "/pfe.plt",
             test_pots * "/ffe.plt",
             test_pots * "/F_fe.plt")

## equilibrate a unit cell
fe1 = bulk(:Fe)
set_calculator!(fe1, eam_Fe)
variablecell!(fe1)
minimise!(fe1)

## get the force constants
fcm = ForceConstantMatrix1(eam_Fe, fe1, h = 1e-5)

println("compute exact hessian on a large cell")
at = fe1 * 12
fixedcell!(at)
set_calculator!(at, eam_Fe)
H = JuLIP.hessian_pos(eam_Fe, at)
Hvec = JuLIP.hessian(at)

println("random virtual displacement => compare H * u and fcm * u")
U = rand(JVecF, length(at))
Uvec = mat(U)[:]
V = vecs(Hvec * Uvec)
println("sanity check")
println(@test V[123] ≈ sum( H[123,n] * U[n] for n = 1:length(U) ))

## now multiply FCM * U
Vfcm = fcm * (at, U)
err_HxU = maximum(norm.(Vfcm - V)) / maximum(norm.(V))
@show err_HxU
println(@test err_HxU < 1e-4)
@warn("""The tolerance for this test has been increased from 1e-7 to 1e-4; see
         also issue #17""")

##

println("Check that individual blocks match")
X = positions(at)
x̄ = sum(X) / length(X)
r = [norm(x-x̄) for x in X]
n0 = findall( r .== minimum(r) )[1]
all_found = true
all_passed = true
for n = 1:length(at)
   n == n0 && continue
   norm(H[n0,n]) < 1e-12 && continue
   R = X[n] - X[n0]
   found = false
   for m = 1:length(fcm.R)
      if norm(R - fcm.R[m]) < 1e-5
         found = true
         if norm(H[n0, n] - fcm.H[m]) < 1e-4
            print("+")
         else
            global all_passed = false
            print("-")
         end
         break
      end
   end
   if !found
      global all_found = false
      print("x")
   end
end
println(@test all_passed)
println(@test all_found )
@warn("several tolerances have been increased for this test, cf #17")


text: include("Type/Type.jl")


text:
"""
```
function interp(x::Vector{T}, x_knot::Vector{T}) where T<:AbstractFloat
```
Compute linear interpolation over 1D grid. It extrapolates outside of
knot points. (SeHyoun Ahn, March 2017)

# Arguments:
- `x_fine`: x grid to interpolate to
- `x_knot`: x knot points to interpolate from

# Output:
`V`: matrix giving interpolation values

# Example:
```
x_fine = linspace(-1,1,300)'
x_knot = linspace(0,2,20)'
V = interp(x_fine, x_knot)
plot(V * (x_knot.^5 + x_knot * 0.5))
```
"""
function interp(x::Vector{T}, x_knot::Vector{T}) where {T<:AbstractFloat}
    loc = sum(broadcast(-, x, x_knot') .>= 0, dims=2)
    loc = max.(min.(loc, length(x_knot) - 1), 1)

    t = (x - x_knot[loc]) ./ (x_knot[loc+1] - x_knot[loc])
    ind_x = 1:length(x)

    i = vec(repeat(ind_x,2,1))
    j = [loc; loc + 1]
    v = [(1 - t); t]

    return SparseArrays.sparse(i, j, v, length(x), length(x_knot))
end

"""
```
function interp(x_fine::Vector{T}, y_fine::Vector{T},
                x_knot::Vector{T}, y_knot::Vector{T}) where {T<:Float64}
```
Compute linear interpolation over 2D grid. It extrapolates outside of
knot points. (SeHyoun Ahn, March 2017)

# Arguments
 - `x_fine::Vector{T}`: x grid to interpolate to
 - `y_fine::Vector{T}`: y grid to interpolate to
 - `x_knot::Vector{T}`: x knot points to interpolate from
 - `y_knot::Vector{T}`: y knot points to interpolate from

# Output
 - `V`: sparse matrix giving interpolation values

# Example
```
x_fine = linspace(0,2,80)'
y_fine = linspace(-1,1,100)'
x_knot = linspace(0,1,10)'
y_knot = linspace(0,1,10)'

V = interp(x_fine, y_fine, x_knot, y_knot)

z = broadcast(+, x_knot.^3 + exp(-x_knot), (y_knot' - 0.5).^2)
surf(x_fine,y_fine,reshape(V * vec(z), 80, 100)')
```
"""
function interp(x_fine::Vector{T}, y_fine::Vector{T},
                x_knot::Vector{T}, y_knot::Vector{T}) where {T<:Float64}

    n_x      = length(x_knot)
    ind_fine = 1:length(x_fine) * length(y_fine)

    loc_x = sum(broadcast(-, x_fine, x_knot') .>= 0, dims=2)
    loc_x = vec(max.(min.(loc_x, n_x - 1), 1))

    loc_y = sum(broadcast(-, y_fine, y_knot') .>= 0, dims=2)
    loc_y = max.(min.(loc_y, length(y_knot) - 1), 1)'

    t_x = (x_fine .- x_knot[loc_x]) ./ (x_knot[loc_x .+ 1] .- x_knot[loc_x])
    t_y = ((y_fine .- vec(y_knot[loc_y])) ./
           (vec(y_knot[loc_y .+ 1]) .- vec(y_knot[loc_y])))'

    j_SW = broadcast(+, loc_x,      n_x.*(loc_y .- 1))
    j_SE = broadcast(+, loc_x .+ 1, n_x.*(loc_y .- 1))
    j_NW = broadcast(+, loc_x,      n_x.*(loc_y))
    j_NE = broadcast(+, loc_x .+ 1, n_x.*(loc_y))
    v_SW = broadcast(*, (1 .- t_x), (1 .- t_y))
    v_SE = broadcast(*, (t_x),      (1 .- t_y))
    v_NW = broadcast(*, (1 .- t_x), (t_y))
    v_NE = broadcast(*, (t_x),      (t_y))

    i = vec(repeat(ind_fine, 4, 1))
    j = [vec(j_SW); vec(j_SE); vec(j_NW); vec(j_NE)]
    v = [vec(v_SW); vec(v_SE); vec(v_NW); vec(v_NE)]

    return SparseArrays.sparse(i, j, v, length(x_fine) * length(y_fine), n_x * length(y_knot))
end


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

@testset "Operating system predicates" begin
    @test !Sys.isunix(:Windows)
    @test !Sys.islinux(:Windows)
    @test Sys.islinux(:Linux)
    @test Sys.iswindows(:Windows)
    @test Sys.iswindows(:NT)
    @test !Sys.iswindows(:Darwin)
    @test Sys.isapple(:Darwin)
    @test Sys.isapple(:Apple)
    @test !Sys.isapple(:Windows)
    @test Sys.isunix(:Darwin)
    @test Sys.isunix(:FreeBSD)
    for bsd in (:FreeBSD, :OpenBSD, :NetBSD, :DragonFly)
        f = Symbol("is", lowercase(String(bsd)))
        q = QuoteNode(bsd)
        @eval begin
            @test Sys.$f($q)
            @test Sys.isbsd($q)
            @test Sys.isunix($q)
            @test !Sys.isapple($q)
        end
    end
    @test_throws ArgumentError Sys.isunix(:BeOS)
    if !Sys.iswindows()
        @test Sys.windows_version() == v"0.0.0"
    else
        @test Sys.windows_version() >= v"1.0.0-"
    end
end

@testset "@static" begin
    @test (@static true ? 1 : 2) === 1
    @test (@static false ? 1 : 2) === 2
    @test (@static if true 1 end) === 1
    @test (@static if false 1 end) === nothing
    @test (@static true && 1) === 1
    @test (@static false && 1) === false
    @test (@static true || 1) === true
    @test (@static false || 1) === 1
    @test (@static if false 1 elseif true 2 end) === 2
    @test (@static if false 1 elseif false 2 end) === nothing
    @test (@static if false 1 elseif false 2 else 3 end) === 3
    @test (@static if false 1 elseif false 2 elseif true && false 3 else 4 end) === 4
    @test (@static if false 1 elseif false 2 elseif true && false 3 end) === nothing
end

if Sys.iswindows()
    @testset "path variables use correct path delimiters on windows" begin
        for path in (Base.SYSCONFDIR, Base.DATAROOTDIR, Base.DOCDIR,
                     Base.LIBDIR, Base.PRIVATE_LIBDIR, Base.INCLUDEDIR)
            @test !occursin("/", path)
            @test !occursin("\\\\", path)
        end
    end
end

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
export Uniform

doc"""
    Uniform(a, b)
    Uniform(b)

Generator of ndarray with a uniform distribution.

# Arguments
* a: Lower bound of the range of random values.
* b: Upper bound of the range of random values.
"""
struct Uniform
    a
    b

    function Uniform(a, b)
        @assert a <= b
        new(a, b)
    end
end

Uniform(b) = Uniform(-b, b)

function (u::Uniform)(::Type{T}, dims::Int...) where T
    r = rand(T, dims)
    r .*= (u.b - u.a)
    r .-= u.a
    r
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1HostAlias <: SwaggerModel
    hostnames::Any # spec type: Union{ Nothing, Vector{String} } # spec name: hostnames
    ip::Any # spec type: Union{ Nothing, String } # spec name: ip

    function IoK8sApiCoreV1HostAlias(;hostnames=nothing, ip=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1HostAlias, Symbol("hostnames"), hostnames)
        setfield!(o, Symbol("hostnames"), hostnames)
        validate_property(IoK8sApiCoreV1HostAlias, Symbol("ip"), ip)
        setfield!(o, Symbol("ip"), ip)
        o
    end
end # type IoK8sApiCoreV1HostAlias

const _property_map_IoK8sApiCoreV1HostAlias = Dict{Symbol,Symbol}(Symbol("hostnames")=>Symbol("hostnames"), Symbol("ip")=>Symbol("ip"))
const _property_types_IoK8sApiCoreV1HostAlias = Dict{Symbol,String}(Symbol("hostnames")=>"Vector{String}", Symbol("ip")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1HostAlias }) = collect(keys(_property_map_IoK8sApiCoreV1HostAlias))
Swagger.property_type(::Type{ IoK8sApiCoreV1HostAlias }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1HostAlias[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1HostAlias }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1HostAlias[property_name]

function check_required(o::IoK8sApiCoreV1HostAlias)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1HostAlias }, name::Symbol, val)
end