text:
using LinearAlgebra
using CVChannel
using Test

"""
This verifies the multiplicativity of the phase rotated amplitude channel family.

The channel family is parameterized by one parameters. We scan over this parameter
and check that it is multiplicative at each point.
...
"""

println("First we define the channel and verify it works.")
function phaseRotAmpDamp(ρ :: Matrix{<:Number}, ε :: Union{Int,Float64}) :: Matrix{<:Number}
    K0 = [1 0 ; 0 sqrt(1-ε)]
    K1 = 1/sqrt(2)*[1 exp(1im*pi/4) ; exp(-1im*pi/4) 1]*[0 sqrt(ε) ; 0 0]

    return K0*ρ*K0' + K1*ρ*K1'
end

@testset "Verify the phase rotated Amplitude Damping Channel Works" begin
    ρ1 = [3/4 0 ; 0 1/4]
    ρ2 = 1/2*[1 -1im ; 1im 1]
    ε1 = 0.3; ε2 = 0.7;

    target1(ε) = [
        3/4+1/8*ε  1/8*exp(1im*pi/4)*ε
        1/8*exp(-1im*pi/4)*ε 1/4-1/8*ε
    ]

    target2(ε) = [
        1/2+1/4*ε -1im/2*sqrt(1-ε)+1/4*exp(1im*pi/4)*ε
        1im/2*sqrt(1-ε)+1/4*exp(-1im*pi/4)*ε 1/2-1/4*ε
    ]

    @test isapprox(phaseRotAmpDamp(ρ1,ε1),target1(ε1), atol=1e-6)
    @test isapprox(phaseRotAmpDamp(ρ1,ε2),target1(ε2), atol=1e-6)
    @test isapprox(phaseRotAmpDamp(ρ2,ε1),target2(ε1), atol=1e-6)
    @test isapprox(phaseRotAmpDamp(ρ2,ε2),target2(ε2), atol=1e-6)
end

println("Now we are going to verify the multiplicativity of the channel")
@testset "Verify Multiplicativity of the Channel" begin
    ε_vals = [0:0.01:1;];
    ctr = 1;
    results = zeros(size(ε_vals)[1],4)
    println("\nBeginning to get the data...")
    println("\nPlease note it starts slow calling the solver.")
    for ε_val in ε_vals
        if (ε_val != 0)&&(Int(floor(ε_val*100)) % 5 == 0)
            println("\nNow evaluating for ε = ", ε_val)
        end
        phaseRotAmpDampVarep(ρ) = phaseRotAmpDamp(ρ,ε_val)
        orig_choi = choi(phaseRotAmpDampVarep,2,2)
        test1 = pptCVDual(orig_choi,2,2)
        kron_par_choi = kron(orig_choi,orig_choi)
        par_choi = permuteSubsystems(kron_par_choi,[1,3,2,4],[2,2,2,2])
        test2 = pptCVPrimal(par_choi,4,4)
        results[ctr,:] = [ε_val, test1[1], test2[1], test2[1]-test1[1]^2]
        ctr = ctr + 1
    end
    println("\ncolumn labels at bottom")
    show(stdout, "text/plain", results)
    println("\n p     cv(N)   cv(N^2)  diff")
    println("\nThus by looking at the diff column, we see what we were verifying")
    println("up to a numerical accuracy of 1e-5 (some computers can do 5e-6).")
    @test all(result -> result > 0 || isapprox(result, 0, atol=1e-5), results[:,4])
end

max_forks_repo_name: ChitambarLab/cv-channel


text: include(joinpath("partials", "item.json.jl"))(@vars(:todo))

ext: jl

max_stars_repo_name: essenciary/genie-todo-mvc

max_issues_repo_name: essenciary/genie-todo-mvc

max_forks_repo_name: essenciary/genie-todo-mvc


text:
# Autogenerated wrapper script for rmapi_jll for aarch64-linux-musl
export rmapi

JLLWrappers.@generate_wrapper_header("rmapi")
JLLWrappers.@declare_executable_product(rmapi)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        rmapi,
        "bin/rmapi",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
Created on Friday 26 July 2019
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Standard matrix factorization algorithms applied on Kronecker systems.
=#

abstract type FactorizedKronecker{T} <: AbstractKroneckerProduct{T} end

# CHOLESKY DECOMPOSITION
# ----------------------

import LinearAlgebra: Cholesky, cholesky

struct CholeskyKronecker{T<:Union{Cholesky,FactorizedKronecker},S<:Union{Cholesky,FactorizedKronecker}} <: FactorizedKronecker{Float64}
    A::T
    B::S
end

issquare(C::Cholesky) = true

function Base.getproperty(C::CholeskyKronecker, d::Symbol)
    if d in [:U, :L, :UL]
        return kronecker(getproperty(C.A, d), getproperty(C.B, d))
    elseif d in [:A, :B]
        return getfield(C, d)
    else
        throw(ArgumentError("Attribute :$d not supported (only :A, :B, :UL, :U or :L)"))
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, C::CholeskyKronecker)
    summary(io, C); println(io)
    println(io, "U factor:")
    show(io, mime, getproperty(C, :U))
end

"""
    cholesky(K::AbstractKroneckerProduct; check = true)

Wrapper around `cholesky` from the `LinearAlgebra` package. Performs Cholesky
on the matrices of a `AbstractKroneckerProduct` instances and returns a
`CholeskyKronecker` type. Similar to `Cholesky`, `size`, `\\`, `inv`, `det`,
and `logdet` are overloaded to efficiently work with this type.
"""
function cholesky(K::AbstractKroneckerProduct; check = true)
    checksquare(K)
    A, B = getmatrices(K)
    return CholeskyKronecker(cholesky(A, check=check),
                            cholesky(B, check=check))
end


text:
#!/usr/bin/env julia
include("../src/Reachability.jl")
using Base.Test, LazySets, Reachability, MathematicalSystems

include("Systems/alltests.jl")
include("ReachSets/alltests.jl")
include("Reachability/alltests.jl")
include("Properties/alltests.jl")


text:
using LinearAlgebra
using StaticArrays
export factorial2, dist2, ipairs, rpairs, spairs, iindex

@inline factorial2(n::Int64) = prod(n:-2:1) # double factorial !!
@inline dist2(dx,dy,dz) = dx*dx+dy*dy+dz*dz # TODO: use hypot()^2?
@inline dist2(dxyz) = dot(dxyz,dxyz)
@inline dist2(xyz1,xyz2) = dist2(xyz1-xyz2)

ipairs(n::Int64) = ((i, j) for i = 1:n for j = 1:i)
rpairs(n::Int64) = ((i,j) for i in 1:n for j in 1:n) # rectangular option to old pairs
spairs(n::Int64) = ((i, j) for i = 1:n for j = 1:(i-1)) # subdiagonal option to old pairs

triangle(i::Int64) = div(i*(i+1),2)
triangle(i::Int64,j::Int64) = i<j ? triangle(j-1)+i : triangle(i-1)+j

@inline iiterator(n::Int64) = ((i,j,k,l) for (i,j) in ipairs(n) for (k,l) in ipairs(n) if triangle(i,j) <= triangle(k,l))

@inline iindex(i::Int64,j::Int64,k::Int64,l::Int64) = triangle(triangle(i,j),triangle(k,l))

struct Val{x} end
Val(x) = Val{x}()


text:
"""
    spectracn(spectra; kwargs...)

Draw a spectra copy number plot for a 2D k-mer frequency spectra.
The spectracn plot consists of a set of differently coloured, stacked barplots.

# Arguments
- `spectra`: A 2D Kmer frequency spectra (`KmerFrequencySpectra{2}`)
"""
@recipe(SpectraCN, spectra) do scene
    default_theme(scene, BarPlot)
end

const CN_COLOUR = [:black, :red, :purple, :green, :blue, :yellow, :orange]

function AbstractPlotting.plot!(plot::SpectraCN)
    # Lower to multiple barplots/plots here.
    spct = to_value(plot[:spectra])
    heights = cumsum(_find_plottable_subset(spct); dims = 2)
    @inbounds for col in size(heights)[2]:-1:1
        barplot!(plot, [Point2f0(i - 1, j) for (i, j) in enumerate(heights[:,col])], color = CN_COLOUR[col])
    end
    # TODO: Figure out how to add a colour legend and perhaps axis labels.
    #xlabel!(plot, "K-mer multiplicity")
    #ylabel!(plot, "Number of distinct k-mers")
    #plot.parent[Axis].attributes[:names][:axisnames] = ("KmerMultiplicity", "Number of distinct k-mers")
    return plot
end


text:
#Write a program that prints a multiplication table for numbers up to 12.

using Printf

println(" X |   1   2   3   4   5   6   7   8   9  10  11  12")
println("---+------------------------------------------------")

for i=1:12, j=0:12
    if j == 0
        @printf("%2d | ", i)
    elseif i <= j
        @printf("%3d%c", i * j, j == 12 ? '\n' : ' ')
    else
        print("    ")
    end
end

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
if !@isdefined USE_GPU
    const USE_GPU = length(ARGS) > 0 && ARGS[1] == "gpu"
end
using ReferenceTests, Test
using CUDA
using LinearAlgebra
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
if !ParallelStencil.is_initialized()
    @static if USE_GPU
        @init_parallel_stencil(CUDA, Float64, 2)
    else
        @init_parallel_stencil(Threads, Float64, 2)
    end
end
include("../scripts-part2/multigrid.jl")

@testset "Test backslash" begin
    nx = 6
    ny = 6
    h = 1 / (ny - 1)
    x = rand((nx-2)*(ny-2))
    A = stencil_5pt(nx-2, ny-2) / h^2
    b = A * x
    xhat = A \ b
    @show norm(x - xhat) / norm(x)
    @test norm(x - xhat) / norm(x) < 1e-10
end

@testset "Test Multigrid with policy=$(execution_policy) and coarse_solver=$(solver) and coarse_solve_size = $((2^l)+1) and nx=ny=$((2^k)+1) on $(USE_GPU ? "GPU" : "CPU")." for execution_policy in [parallel, parallel_shmem], l in 2:3, k in 7:10, solver in [jacobi, conjugate_gradient]
    n = (2^k)+1
    h = 1 / (n - 1)
    c = 0.0
    tol = 1e-6

    opt = MGOpt()
    opt.execution_policy = execution_policy
    opt.coarse_solve_size = (2^l)+1
    opt.coarse_solver = solver

    inn = CartesianIndices((2:n-1, 2:n-1))

    xref = zeros(n, n)
    xref[inn] .= rand(n-2, n-2)
    b = zeros(n, n)
    b_ = zeros(n-2, n-2)
    A = stencil_5pt(n-2, n-2) / h^2
    b_[:] .= A * xref[inn][:]
    b[inn] .= b_

    xhat = @zeros(n, n)
    r_rms = MGsolve_2DPoisson!(xhat, Data.Array(b), h, c, tol, 20, false; opt=opt, verbose=true)
    @synchronize()
    xhat = Array(xhat)

    @test r_rms < tol * sqrt( sum(b.^2) / (n*n) )
    #@test norm(xref - xhat) / norm(xref) < tol
end

@testset "Test Jacobi solver" for execution_policy in [parallel, parallel_shmem]
    n = 33
    h = 1 / (n - 1)
    c = 0.0
    alpha = 1.0
    tol = 1e-6
    Nmax = 10000

    inn = CartesianIndices((2:n-1, 2:n-1))

    # CPU
    xref = rand(n, n)
    xref[1, :] .= 0.
    xref[n, :] .= 0.
    xref[:, 1] .= 0.
    xref[:, n] .= 0.

    b = zeros(n, n)
    b_ = zeros(n-2, n-2)

    A = stencil_5pt(n-2, n-2) / h^2
    b_[:] .= A * xref[inn][:]
    b[inn] .= b_

    # XPU
    res_buf = @zeros(n, n)
    b = Data.Array(b)
    xhat = @zeros(n, n)
    tolb = tol * sqrt(sum(b.^2)/(n*n))
    for i = 1 : Nmax
        res_rms = iteration_2DPoisson!(xhat, b, h, c, res_buf, execution_policy)
        @synchronize
        if res_rms < tolb
            println("iter = $(i)")
            break
        end
    end
    xhat = Array(xhat)

    @test norm(xref - xhat) / norm(xref) < tolb
end

@testset "Test residual_2DPoisson" for execution_policy in [parallel, parallel_shmem]
    # CPU matrix
    # XPU residual
    n = 64
    h = 1 / (n - 1)
    c = 3.1415
    inn = CartesianIndices((2:n-1, 2:n-1))

    # XPU
    u = @rand(n, n)
    u[1, :] .= 0.
    u[n, :] .= 0.
    u[:, 1] .= 0.
    u[:, n] .= 0.
    f = @rand(n, n)
    res = @zeros(n, n)

    residual_2DPoisson_wrapper!(u, f, h, c, res, execution_policy)
    @synchronize()

    res_cpu = Array(res)

    # CPU

    u_cpu = Array(u)
    u_ = u_cpu[inn]

    f_cpu = Array(f)
    f_ = f_cpu[inn]

    A = stencil_5pt(n-2, n-2) / h^2 - c * I

    res_ = zeros(n-2, n-2)
    res_[:] .= A * u_[:] - f_[:]

    @test res_ ≈ res_cpu[inn]
end

@ParallelStencil.reset_parallel_stencil()


text:
using MagmaThermoKinematics
using MagmaThermoKinematics.Diffusion2D
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
#using Plots
using CairoMakie
using GeoParams
using Printf

# Initialize
@init_parallel_stencil(Threads, Float64, 2);    # initialize parallel stencil in 2D

#------------------------------------------------------------------------------------------
@views function MainCode_2D();

K = u"K"    # unsure why this is not available from GeoParams

CharDim      = GEO_units(length=1000km, temperature=1000C, stress=10MPa, viscosity=1e20Pas);

# Define material parameters for the simulation. By using GeoParams
MatParam     = Array{MaterialParams, 1}(undef, 2);
MatParam[1]  = SetMaterialParams(Name="Host Rock", Phase=1,
                                 Density    = ConstantDensity(ρ=2700/m^3),
                          EnergySourceTerms = ConstantLatentHeat(Q_L=313/kg),
                               #Conductivity = ConstantConductivity(k=1.5Watt/K/m),     # in case we use constant k
                               Conductivity = T_Conductivity_Whittacker(),
                               HeatCapacity = ConstantHeatCapacity(cp=1000/kg/K),
                                    Melting = MeltingParam_Caricchi())

MatParam[2]  = SetMaterialParams(Name="Melt", Phase=2,
                                    Density    = ConstantDensity(ρ=2700kg/m^3),
                             EnergySourceTerms = ConstantLatentHeat(Q_L=313kJ/kg),
                                  Conductivity = ConstantConductivity(k=1.2Watt/K/m),
                                  HeatCapacity = ConstantHeatCapacity(cp=1000/kg/K),
                                       Melting = MeltingParam_Caricchi())

# These
SecYear = 3600*24*365.25
kyr = 1000*SecYear
km³ = 1000^3

# Model parameters
W,H                     =   25e3, 20e3;                 # Width, Height
ρ                       =   2700;                       # Density
cp                      =   1000;                       # Heat capacity
k_rock, k_magma         =   1.5, 1.2;                   # Thermal conductivity of host rock & magma
La                      =   350e3;                      # Latent heat J/kg/K
GeoT                    =   50.0/1000;                  # Geothermal gradient [K/km]
W_in, H_in              =   20e3,    100;               # Width and thickness of each dike that is intruded

T_in                    =   1115;                       # Intrusion temperature
maxTime                 =   1250kyr;                    # Maximum simulation time in kyrs

TotalVolume             =   1580e3^3
maxThickness            =   7000                        # max thickness of dike
nInjections             =   maxThickness/H_in
InjectionInterval       =   maxTime/nInjections;        # Inject a new dike every X kyrs

DikeType                =   "SquareDike"              # Type to be injected ("ElasticDike","SquareDike")
#DikeType                =   "SquareDike_TopAccretion"   # Type to be injected ("ElasticDike","SquareDike", "SquareDike_TopAccretion")
#DikeType                =   "ElasticDike"   # Type to be injected ("ElasticDike","SquareDike", "SquareDike_TopAccretion")

cen                     =   [0.; -6.5e3 - 100/2];       # Center of dike

Nx, Nz                  =   200, 250;               # resolution (grid cells in x,z direction)
dx,dz                   =   W/(Nx-1), H/(Nz-1);     # grid size
κ                       =   k_rock./(ρ*cp);         # thermal diffusivity
dt                      =   min(dx^2, dz^2)./κ/20;  # stable timestep (required for explicit FD)
nt::Int64               =   floor(maxTime/dt);      # number of required timesteps
nTr_dike                =   300;                    # number of tracers inserted per dike

@show nt maxTime dt

# Array initializations
T                       =   @zeros(Nx,   Nz);
Kc                      =   @ones(Nx,    Nz)*k_rock;
Rho                     =   @ones(Nx,    Nz)*ρ;
Cp                      =   @ones(Nx,    Nz)*cp;

# Work array initialization
Tnew, qr,qz, Kr, Kz     =   @zeros(Nx,   Nz), @zeros(Nx-1, Nz),     @zeros(Nx,   Nz-1), @zeros(Nx-1, Nz), @zeros(Nx,   Nz-1)    # thermal solver
R,Rc,Z                  =   @zeros(Nx,   Nz), @zeros(Nx-1, Nz-1),   @zeros(Nx,   Nz)    # 2D gridpoints
Phi_o, Phi, dPhi_dt     =   @zeros(Nx,   Nz), @zeros(Nx,   Nz  ),   @zeros(Nx,   Nz)    # solid fraction

# Set up model geometry & initial T structure
x,z                     =   (0:Nx-1)*dx, (-(Nz-1):0)*dz;                            # 1-D coordinate arrays
crd                     =   collect(Iterators.product(x,z))                         # Generate coordinates from 1D coordinate vectors
R,Z                     =   (x->x[1]).(crd), (x->x[2]).(crd);                       # Transfer coords to 3D arrays
Rc                      =   (R[2:end,:] + R[1:end-1,:])/2
Grid                    =   (x,z);                                                  # Grid
Tracers                 =   StructArray{Tracer}(undef, 1)                           # Initialize tracers
dike                    =   Dike(W=W_in,H=H_in,Type=DikeType,T=T_in);               # "Reference" dike with given thickness,radius and T
T                       .=   -Z.*GeoT;                                              # Initial (linear) temperature profile
Phi, dPhi_dt            =   SolidFraction(T, Phi_o, dt);                            # Compute solid fraction

# Preparation of visualisation
#ENV["GKSwstype"]="nul"; if isdir("viz2D_out")==false mkdir("viz2D_out") end; loadpath = "./viz2D_out/"; anim = Animation(loadpath,String[])

# Create main figure, which will be updated every time the obsevable changes
fig = Figure(resolution = (960,450))
Phi_melt    =   1.0 .- Phi;

time=0.0
ax1 = Axis(fig, aspect = 1, xlabel="width [km]", ylabel="Depth [km]", title=" $(round(time/kyr*100)/100) kyrs ")
ax2 = Axis(fig, aspect = 1, xlabel="width [km]")
cmap1 = :rainbow1
#cmap1 = :vik   # better because perceptionally uniform
cmap2 = :acton

pltobj1 = heatmap!(ax1, x/1e3, z/1e3, T, colorrange = (0,1000), colormap = cmap1)
contour!(ax1, x/1e3, z/1e3, T,levels = 125:125:1000,color = :black,linewidth = 0.85)
pltobj2 = heatmap!(ax2, x/1e3, z/1e3, Phi_melt, colorrange = (0,1),colormap = cmap2)

cbar1 = Colorbar(fig, pltobj1, label="Temperature", ticks=125:125:1000)
cbar2 = Colorbar(fig, pltobj2, label="Melt fraction")

o = Observable(1)
on(o) do o
    Phi_melt    =   1.0 .- Phi;

    ax1.title=" $(round(time/kyr*100)/100) kyrs "   # update title
    pltobj1 = heatmap!(ax1, x/1e3, z/1e3, T, colorrange = (0,1000), colormap = cmap1)
              contour!(ax1, x/1e3, z/1e3, T,levels = 100:100:1000,color = :black,linewidth = 0.85)
    pltobj2 = heatmap!(ax2, x/1e3, z/1e3, Phi_melt, colorrange = (0,1),colormap = cmap2)
    fig[1,1] = ax1
    fig[1,2] = cbar1
    fig[1,3] = ax2
    fig[1,4] = cbar2
    display(fig)

    save("./viz2D_out/$(num+10000).png", fig)    # also save as separeta
end

num=0;
time, dike_inj, InjectVol, Time_vec,Melt_Time = 0.0, 0.0, 0.0,zeros(nt,1),zeros(nt,1);

record(fig, "sill_intrusion.gif") do io  # we will save the animation in the gif

    for it = 1:nt   # Time loop

        if floor(time/InjectionInterval)> dike_inj              # Add new dike every X years
            dike_inj   =   floor(time/InjectionInterval)        # Keeps track on what was injected already
            Angle_rand =   0;
            dike       =   Dike(dike, Center=cen[:],Angle=[Angle_rand]);                # Specify dike with random location/angle but fixed size/T
            Tracers, T, Vol     =   InjectDike(Tracers, T, Grid, dike, nTr_dike);       # Add dike, move hostrocks
            InjectVol           +=  Vol                                                 # Keep track of injected volume
            Qrate           = InjectVol/time
            Qrate_km3_yr    = Qrate*SecYear/(1e3^3)
            Qrate_km3_yr_km2    = Qrate_km3_yr/(pi*(W_in/2/1e3)^2)

            @printf "Added new dike; total injected magma volume = %.2f km³; rate Q=%.2f m³s⁻¹ = %.2e km³yr⁻¹ = %.2e km³yr⁻¹km⁻² \n" InjectVol/km³ Qrate Qrate_km3_yr Qrate_km3_yr_km2

        #   println("Added new dike; total injected magma volume = $(round(InjectVol/km³,digits=2)) km³; rate Q=$(round(Qrate,digits=2)) m³/s, Q=$(Qrate_km3_yr) km³/yr = $(Qrate_km3_yr_km2)km³/yr/km^2")
        end

        T_K         =  (T .+ 273.15).*K
        Phi_melt    =   ComputeMeltingParam(0,T_K,MatParam[1].Melting[1])
        Phi         =   1.0 .- Phi_melt;
        dPhi_dt     =   (Phi-Phi_o)./dt
    #  Phi, dPhi_dt     =  SolidFraction(T, Phi_o, dt);                                # Compute solid fraction

        Kc          =   ComputeConductivity(0,T_K,MatParam[1].Conductivity[1])           # Thermal conductivity
        Kc          =   ustrip.(Kc)

        # Perform a diffusion step
        #@parallel diffusion2D_step!(Tnew, T, qx, qz, K, Kr, Kz, Rho, Cp, dt, dr, dz, La, dPhi_dt);
        @parallel diffusion2D_AxiSymm_step!(Tnew, T, R, Rc, qr, qz, Kc, Kr, Kz, Rho, Cp, dt, dx, dz, La, dPhi_dt)

        # set lateral boundary conditions (flux-free)
        @parallel (1:size(T,2)) bc2D_x!(Tnew);

        #Tnew[:,1] .= GeoT*H;                               # fixed bottom temperature (constant)
        @parallel (1:size(T,1)) bc2D_z_bottom!(Tnew);       # flux-free bottom BC

        # Fixed top BC
        Tnew[:,end] .= 0.0;                                 # fixed top temperature (constant)

        Tracers             =   UpdateTracers(Tracers, Grid, Tnew, Phi);                            # Update info on tracers
        T, Tnew             =   Tnew, T;                                                            # Update temperature
        time                =   time + dt;                                                          # Keep track of evolved time
        Melt_Time[it]       =   sum( 1.0 .- Phi)/(Nx*Nz)                                            # Melt fraction in crust
        Time_vec[it]        =   time;                                                               # Vector with time

        if mod(it,1000)==0  # Visualisation
            num = num+1;

            o[] = num;      # Creates the updated figure
            recordframe!(io)
            println(" Timestep $it = $(round(time/kyr*100)/100) kyrs, maxT = $(maximum(T)) C")
        end
    end
end
#save("test.gif", stream)

#gif(anim, "Example2D.gif", fps = 15)   # create gif animation
return Time_vec, Melt_Time;
end # end of main function

Time_vec,Melt_Time = MainCode_2D(); # start the main code
#plot(Time_vec/kyr, Melt_Time, xlabel="Time [kyrs]", ylabel="Fraction of crust that is molten", label=:none); png("Time_vs_Melt_Example2D") #Create plot


text:
using Compat
using VoronoiDelaunay
using Images

include("triang.jl")
include("polyline.jl")
include("prepared.jl")
include("simple.jl")
include("view.jl")


text:
Rt2T(R,t) = [R t; 0 0 0 1]
T2R(T::AbstractMatrix) = T[1:3,1:3]
T2t(T::AbstractMatrix) = T[1:3,4]
skewcoords(R) = [R[3,2];R[1,3];R[2,1]]
twistcoords(xi) = [xi[1:3, 4]; skewcoords(xi[1:3, 1:3])]
@inline skew(s)::AbstractMatrix{eltype(s)} = [0 -s[3] s[2];s[3] 0 -s[1]; -s[2] s[1] 0]
@inline function skew!(R::T,s)::T where T
    R[1,1] = 0
    R[2,1] = s[3]
    R[3,1] = -s[2]
    R[1,2] = -s[3]
    R[2,2] = 0
    R[3,2] = s[1]
    R[1,3] = s[2]
    R[2,3] = -s[1]
    R[3,3] = 0
    return R
end
skew4(s)::AbstractMatrix{eltype(s)} = [skew(s[4:6]) s[1:3]; 0 0 0 0]

# expω(w,q=1) = I + sin(norm(w)*q)/norm(w)*skew(w) + (1-cos(norm(w)*q))/norm(w)^2*skew(w)^2 # verified to work

function expω(w,q=1)
    nw = norm(w)
    if nw < 1e-12 # TODO: Use TaylorSeries.jl to approximate this for small nw, I verified TaylorSeries.jl to work well for this
        I + q*skew(w)
    else
        I + sin(nw*q)/nw*skew(w) + (1-cos(nw*q))/nw^2*skew(w)^2
    end
end

function expξ2(xi,q=1.0) # works well for jacobian-ikine
    v = xi[1:3]
    w = xi[4:6]
    nw = norm(w)
    what = skew(w)
    if nw < 1e-12
        eR = I + q*what
        A = q*I
    else
        eR = I + sin(nw*q)/nw*what + (1-cos(nw*q))/nw^2*what^2
        A = q*I + ((1 - cos(nw*q))/(nw^2))*what + ((nw*q - sin(nw*q))/(nw^3))*what*what
    end
    [eR A*v; 0 0 0 1]
end

"""Optimized routine to modify T in place"""
function expξ2!(T,xi,q=1.0)
    w = skew(xi[4:6])
    w2 = w^2
    nw = norm(xi[4:6])
    nw2 = nw^2
    nwq = nw*q
    snwq = sin(nwq)
    cnwq = 1-cos(nwq)
    T[1:3,1:3] = I + snwq/nw*w + cnwq/nw2*w2
    A = q*I + (cnwq/(nw2))*w + ((nwq - snwq)/(nw*nw2))*w2
    T[1:3,4] = A*xi[1:3]
end

"""
`expξ(xi,q=1.0)` calculates the exponential map of a twist with twistcoordinates xi and joint angle q
If no angle is given, q=1 is assumed
"""
function expξ(xi,q=1.0)
    exp(skew4(xi)*q)
end

"""
Calculates the matrix logarithm of a transformation matrix ∈ SE(3)
Does not seem to be very reliable for very small rotations, use logm instead, which is a bit slower.
"""
function logT(T) # Verified to work in the vicinity of I, but not for very small errors
    R = T[1:3,1:3]
    t = T[1:3,4]
    ω = logR(R)
    nω = norm(skewcoords(ω))
    if nω < 1e-3
        fact = 0.0833333
    else
        fact = (2sin(nω)-nω*(1+cos(nω)))/(2*nω^2*sin(nω))
    end
    Ai = I-0.5ω+fact*ω^2
    return [ω Ai*t;0 0 0 0]
end

"""Calculates the matrix logarithm of a rotation matrix ∈ SO(3)"""
function logR(R)
    @assert isrot(R)
    phi = acos((min(tr(R),3)-1)/2)
    if abs(phi) < 1e-10
        ω = (R-R')/2
    else
        ω = phi/(2sin(phi))*(R-R')
    end
    return ω
end

"""Calculates the adjoint of a transformation matrix"""
function ad(T)
    # Computes the adjoint of transformation matrix T
    R = T[1:3,1:3]
    t = skew(T[1:3,4])
    Z = zeros(3,3)
    A = [R t*R; Z R]
    return A
end

"""Calculates the adjoint of the inverse of a tranformation matrix T, which is also the inverse of the adjoint of T"""
function adi(T)
    # Computes the adjoint of transformation matrix inv(T)
    R = T[1:3,1:3]'
    t = skew(-R*T[1:3,4])
    Z = zeros(3,3)
    A = [R t*R; Z R]
    return A
end

"""Inverts a transformation matrix ∈ SE(3)"""
trinv(T) = [T[1:3,1:3]' -T[1:3,1:3]'*T[1:3,4];0 0 0 1]

isrot(R) = det(R[1:3,1:3]) ≈ 1 && norm(R[1:3,1:3]'R[1:3,1:3]-I) < 1e-10
isse3(T) = isrot(T) && T[4,1:4] == [0 0 0 1]

"""`Rangle(R1,R2 = I3,deg = false)` calculates the angle between two rotation matrices"""
function Rangle(R1::AbstractMatrix,R2 = I,deg = false)
    @views cosθ = (tr(R1[1:3,1:3]' * R2[1:3,1:3])-1)/2
    if cosθ > 1
        cosθ = 1
    elseif cosθ < -1
        cosθ = -1
    end
    θ = acos(cosθ)
    if deg
        θ *= 180/pi
    end
    return θ
end

function Rangle(R1::AbstractArray{T,3},R2i = I,deg = false) where T
    N = size(R1,3)
    θ = zeros(N)
    R2 = R2i == I ? I3 : R2i
    for i = 1:N
        cosθ = (tr(R1[1:3,1:3,i]' * R2[1:3,1:3,(size(R2,3)==1 ? 1 : i)])-1)/2
        if cosθ > 1
            cosθ = 1
        elseif cosθ < -1
            cosθ = -1
        end
        θ[i] = acos(cosθ)
        if deg
            θ[i] *= 180/pi
        end
    end
    return θ

end

"""This is a helper method for calibPOE"""
function Ai(q,xi)
    n = size(q,1)
    A = zeros(6,6(n+1))
    qext = [q;1]
    pa = Matrix{Float64}(i, 6, 6)
    for i = 1:n+1
        ω = xi[4:6,i]
        Ω = [skew(xi[4:6,i]) skew(xi[1:3,i]);
        zeros(3,3) skew(xi[4:6,i])] # This might be wrong
        θ = norm(ω)*qext[i]

        Ai = qext[i]*I + (4-θ*sin(θ) - 4cos(θ))/(2norm(ω)^2)* Ω + (4θ-5sin(θ) + θ*cos(θ))/(2norm(ω)^3) *Ω^2 + (2-θ*sin(θ) - 2cos(θ))/(2norm(ω)^4) *Ω^3 + (2θ-3sin(θ) + θ*cos(θ))/(2norm(ω)^5)* Ω^4
        A[:,(i-1)*6+1:6i] = pa*Ai
        pa = pa*ad(expξ(xi[:,i],qext[i])) # Correct order?
    end

    return A
end

"""This is the other helper method for calibPOE"""
function xii(q,xi)
    n = size(q,1)
    A = zeros(6,n)
    pa = Matrix{Float64}(i, 6, 6)
    for i = 1:n
        A[:,i] = pa*xi[:,i]
        pa = pa*ad(expξ(xi[:,i],q[i])) # Correct order?
    end

    return A
end

function prodad(xi,q,n)
    pa = Matrix{Float64}(i, 6, 6)
    for i = 1:n
        pa = pa*ad(expξ(skew4(xi[:,i]),q[i]))
    end
    return pa
end

function conformize(xi)
    xi[4:6] ./= norm(xi[4:6])
    xi[1:3] = xi[1:3] - (xi[4:6]'*xi[1:3])[1]*xi[4:6]
    return xi
end

"""
Calculate the (signed) angle around z axis for two rotation (transformation) matrices.
"""
function xyθ(x1,x2)
    # # Check the angle between the y-axis vectors, should go equally well with x-axis vectors, assuming angles are small
    # cosθ = x1[1:3,2]⋅x2[1:3,2]
    # # The angle error does not have sign after acos. Therefore check the angle between x and y axes and determine the sign from that
    # cosθxy = x1[1:3,1]⋅x2[1:3,2]
    # if cosθ > 1 # To avoid complex results from acos if there is some small numerical error
    #     cosθ = 1
    # end
    # θ = acos(cosθ)
    # if cosθxy < 0 # TODO: this is not properly thought through
    #     θ = -θ
    # end
    # return θ
    R2rpy(x1[1:3,1:3]'x2[1:3,1:3],conv="zyx")[3]
end

"""Takes the DH-parameters or a set of nominal transformation matrices and outputs the joint twists in base frame"""
function DH2twistsPOE(Tn)
    #This is a bit tricky, the last two twists corresponds to [tool T(0)], use with fkinePOE
    n = size(Tn,3)-1
    xi = zeros(6,n+2)
    P = skew4([0,0,0,0,0,1])
    T = I4
    M = I4
    for i = 1:n+1
        Xi = M*P*trinv(M)
        xi[:,i] = twistcoords(Xi)
        T = T*Tn[:,:,i]
        M = M*Tn[:,:,i]
    end
    #xi[:,end-1] = twistcoords(logT(Tn[:,:,end]))
    xi[:,end] = twistcoords(logT(T))
    return xi
end

"""Takes the DH-parameters or a set of nominal transformation matrices and outputs the joint twists in local link frames"""
function DH2twistsLPOE(Tn)
    n = size(Tn,3)
    xi = zeros(6,n)
    P = skew4([0,0,0,0,0,1])
    for i = 1:n
        Xi = trinv(Tn[:,:,i])*P*(Tn[:,:,i])
        xi[:,i] = twistcoords(Xi)
    end
    return xi
end

DH2twistsLPOE(dh::DH) = DH2twistsLPOE(dh2Tn(dh))
DH2twistsPOE(dh::DH) = DH2twistsPOE(dh2Tn(dh))

"""Takes a matrix R ∈ SO(3) or T ∈ SE(3) and makes the rotational part orthonormal"""
function toOrthoNormal!(M)
    R = M[1:3,1:3]
    U,S,V = svd(R)
    a = sign(det(U*V'))
    S = diagm(0=>[1,1,a])
    R = U*S*V'
    M[1:3,1:3] = R
    M
end

function toOrthoNormal(Mi)
    M = deepcopy(Mi)
    R = M[1:3,1:3]
    U,S,V = svd(R)
    a = sign(det(U*V'))
    S = diagm(0=>[1,1,a])
    R = U*S*V'
    M[1:3,1:3] = R
    M
end

"""
`rpy2R(r,p,y,conv="zyx")`
For rpy from ABB robot, use zyx
For rpy from Nikon, use xyz
"""
function rpy2R(r,p,y,conv="zyx")
    # ca = cos(r)
    # cb = cos(p)
    # cg = cos(y)
    # sa = sin(r)
    # sb = sin(p)
    # sg = sin(y)
    # R = [ca*cb*cg-sa*sg    -ca*cb*sg-sa*cg  ca*sb;
    # sa*cb*cg-sa*sg    -sa*cb*sg+ca*cg  sa*sb;
    # -sb*cg            sb*sg            cb]
    if conv == "zyz"
        R = rotz(r)*roty(p)*rotz(y)
    elseif conv == "xyz"
        R = rotx(r)*roty(p)*rotz(y)
    elseif conv == "zyx"
        R = rotz(r)*roty(p)*rotx(y)
    end
end

rpy2R(r,conv="zyx") = rpy2R(r...,conv)

function rotx(t, deg=false)
    if deg
        t = t *pi/180
    end
    ct = cos(t)
    st = sin(t)
    R = [
    1   0    0
    0   ct  -st
    0   st   ct
    ]
end

function roty(t, deg=false)
    if deg
        t = t *pi/180
    end
    ct = cos(t)
    st = sin(t)
    R = [
    ct  0   st
    0   1   0
    -st  0   ct
    ]
end

function rotz(t, deg=false)
    if deg
        t = t *pi/180
    end
    ct = cos(t)
    st = sin(t)
    R = [
    ct  -st  0
    st   ct  0
    0    0   1
    ]
end

"""
`R2rpy(R; conv="xyz", deg = false)`\n
If `conv` is not `xyz`, it will be `zyx`\n
returns a vector ∈ R3 or a matrix ∈ R3×N depending on the dimension of the input
"""
function R2rpy(m::AbstractArray{T,3}; conv="xyz", deg = false) where T
    N = size(m,3)
    rpy = Array(Float64,3,N)
    for i = 1:N
        rpy[:,i] = R2rpy(m[:,:,i], conv=conv, deg=deg)
    end
    return rpy
end

function R2rpy(m::AbstractMatrix; conv="xyz", deg = false)

    rpy = zeros(3)

    if conv == "xyz"
        # XYZ order
        if abs(m[3,3]) < eps() && abs(m[2,3]) < eps()
            # singularity
            rpy[1] = 0;  # roll is zero
            rpy[2] = atan(m[1,3], m[3,3])  # pitch
            rpy[3] = atan(m[2,1], m[2,2])  # yaw is sum of roll+yaw
        else
            rpy[1] = atan(-m[2,3], m[3,3])        # roll
            # compute sin/cos of roll angle
            sr = sin(rpy[1])
            cr = cos(rpy[1])
            rpy[2] = atan(m[1,3], cr * m[3,3] - sr * m[2,3])  # pitch
            rpy[3] = atan(-m[1,2], m[1,1])        # yaw
        end
    else
        # old ZYX order (as per Paul book)
        if abs(m[1,1]) < eps() && abs(m[2,1]) < eps()
            # singularity
            rpy[1] = 0     # roll is zero
            rpy[2] = atan(-m[3,1], m[1,1])  # pitch
            rpy[3] = atan(-m[2,3], m[2,2])  # yaw is difference yaw-roll
        else
            rpy[1] = atan(m[2,1], m[1,1])
            sp = sin(rpy[1])
            cp = cos(rpy[1])
            rpy[2] = atan(-m[3,1], cp * m[1,1] + sp * m[2,1])
            rpy[3] = atan(sp * m[1,3] - cp * m[2,3], cp*m[2,2] - sp*m[1,2])
        end
    end
    if deg
        rpy *= 180/π
    end
    return rpy
end

using Quaternions
import Quaternions.Quaternion
function Quaternion(t::AbstractMatrix{P}) where P
    qs = sqrt(t[1,1]+t[2,2]+t[3,3]+1)/2.0
    kx = t[3,2] - t[2,3]   # Oz - Ay
    ky = t[1,3] - t[3,1]   # Ax - Nz
    kz = t[2,1] - t[1,2]   # Ny - Ox

    if t[1,1] >= t[2,2] && t[1,1] >= t[3,3]
        kx1 = t[1,1] - t[2,2] - t[3,3] + 1 # Nx - Oy - Az + 1
        ky1 = t[2,1] + t[1,2]          # Ny + Ox
        kz1 = t[3,1] + t[1,3]          # Nz + Ax
        add = kx >= 0
    elseif t[2,2] >= t[3,3]
        kx1 = t[2,1] + t[1,2]          # Ny + Ox
        ky1 = t[2,2] - t[1,1] - t[3,3] + 1 # Oy - Nx - Az + 1
        kz1 = t[3,2] + t[2,3]          # Oz + Ay
        add = ky >= 0
    else
        kx1 = t[3,1] + t[1,3]          # Nz + Ax
        ky1 = t[3,2] + t[2,3]          # Oz + Ay
        kz1 = t[3,3] - t[1,1] - t[2,2] + 1 # Az - Nx - Oy + 1
        add = kz >= 0
    end

    if add
        kx = kx + kx1
        ky = ky + ky1
        kz = kz + kz1
    else
        kx = kx - kx1
        ky = ky - ky1
        kz = kz - kz1
    end
    nm = norm([kx, ky, kz])
    if nm == 0
        q = Quaternion(one(P), 0, 0, 0)
    else
        s  = sqrt(1 - qs^2) / nm
        qv = s*[kx, ky, kz]
        q  = Quaternion(qs,qv)
    end
end

function traj2quat(T)
    N = size(T,3)
    Q = Array(Float64,4,N)
    for i = 1:N
        q = Quaternion(T[:,:,i])
        Q[1,i] = q.s
        Q[2,i] = q.v1
        Q[3,i] = q.v2
        Q[4,i] = q.v3
    end
    return Q
end

function centraldiff(v)
    c = size(v,2)
    a1 = [zeros(1,c);diff(v)]
    a2 = [diff(v);zeros(1,c)]
    a = (a1+a2)/2
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI

pathToFMU = joinpath(dirname(@__FILE__), "../model/Dymola/2020x/SpringPendulum1D.fmu")

myFMU = fmiLoad(pathToFMU)

#create an instance and simulate it
comp1 = fmiInstantiate!(myFMU; loggingOn=true)
fmiSetupExperiment(comp1, 0.0)
fmiEnterInitializationMode(comp1)
fmiExitInitializationMode(comp1)

t_start = 0.0
t_stop = 8.0

data1 = fmiSimulateCS(comp1, t_start, t_stop; recordValues=["mass.s"])
fmiPlot(data1)

#create another instance, change the spring stiffness and simulate it
comp2 = fmiInstantiate!(myFMU; loggingOn=true)
fmiSetupExperiment(comp2, 0.0)
fmiEnterInitializationMode(comp2)
springConstant = fmiGetReal(comp2, "spring.c") * 0.1
fmiSetReal(comp2, "spring.c", springConstant)
fmiExitInitializationMode(comp2)
data2 = fmiSimulateCS(comp2, t_start, t_stop; recordValues=["mass.s"])
fmiPlot(data2)

fmiUnload(myFMU)


text:
mutable struct RSWM{T}
  discard_length::T
  adaptivealg::Symbol
end

Base.@pure function RSWM(;
     discard_length=1e-15,
     adaptivealg::Symbol=:RSwM3)
     RSWM{typeof(discard_length)}(discard_length,adaptivealg)
end

adaptive_alg(rswm::RSWM{T}) where {T} = rswm.adaptivealg


text:
function __solve end
function __init end
function solve! end

NO_TSPAN_PROBS = Union{AbstractLinearProblem, AbstractNonlinearProblem,
                       AbstractQuadratureProblem,
                       AbstractSteadyStateProblem,AbstractJumpProblem}

function init_call(_prob,args...;kwargs...)
  if :kwargs ∈ propertynames(_prob)
    __init(_prob,args...;_prob.kwargs...,kwargs...)
  else
    __init(_prob,args...;kwargs...)
  end
end

function init(prob::DEProblem,args...;kwargs...)
  _prob = get_concrete_problem(prob,kwargs)
  if haskey(kwargs,:alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    isadaptive(alg) &&
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    init_call(_prob,alg,args...;kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    isadaptive(alg) &&
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    init_call(_prob,args...;kwargs...)
  else
    init_call(_prob,args...;kwargs...)
  end
end

function solve_call(_prob,args...;merge_callbacks = true, kwargs...)
  if :kwargs ∈ propertynames(_prob)
    if merge_callbacks && haskey(_prob.kwargs,:callback) && haskey(kwargs, :callback)
      kwargs_temp = NamedTuple{Base.diff_names(Base._nt_names(
      values(kwargs)), (:callback,))}(values(kwargs))
      callbacks = NamedTuple{(:callback,)}( [DiffEqBase.CallbackSet(_prob.kwargs[:callback], values(kwargs).callback )] )
      kwargs = merge(kwargs_temp, callbacks)
    end
    kwargs = merge(values(_prob.kwargs), kwargs)
  end

  logger = get(kwargs, :progress, false) ? default_logger(Logging.current_logger()) : nothing
  maybe_with_logger(logger) do
    __solve(_prob,args...; kwargs...)
  end
end

function solve(prob::DEProblem,args...;kwargs...)
  _prob = get_concrete_problem(prob,kwargs)
  if haskey(kwargs,:alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    isadaptive(alg) &&
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    solve_call(_prob,alg,args...;kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    isadaptive(alg) &&
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    solve_call(_prob,args...;kwargs...)
  elseif isempty(args) # Default algorithm handling
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    solve_call(_prob,args...;kwargs...)
  else
    solve_call(_prob,args...;kwargs...)
  end
end

function solve(prob::EnsembleProblem,args...;kwargs...)
  if isempty(args)
    __solve(prob,nothing,args...;kwargs...)
  else
    __solve(prob,args...;kwargs...)
  end
end

function solve(prob::AbstractNoiseProblem,args...;kwargs...)
  __solve(prob,args...;kwargs...)
end

function get_concrete_problem(prob::AbstractJumpProblem,kwargs)
  prob
end

function get_concrete_problem(prob::AbstractSteadyStateProblem, kwargs)
  u0 = get_concrete_u0(prob, Inf, kwargs)
  u0 = promote_u0(u0, prob.p, nothing)
  remake(prob; u0 = u0)
end

function get_concrete_problem(prob::AbstractEnsembleProblem, kwargs)
  prob
end

function DiffEqBase.solve(prob::PDEProblem,alg::DiffEqBase.DEAlgorithm,args...;
                                          kwargs...)
    solve(prob.prob,alg,args...;kwargs...)
end

function discretize end

function get_concrete_problem(prob, kwargs)
  tspan = get_concrete_tspan(prob, kwargs)
  u0 = get_concrete_u0(prob, tspan[1], kwargs)
  u0 = promote_u0(u0, prob.p, tspan[1])
  tspan = promote_tspan(u0, prob.p, tspan, prob, kwargs)
  remake(prob; u0 = u0, tspan = tspan)
end

function get_concrete_problem(prob::DDEProblem, kwargs)
  tspan = get_concrete_tspan(prob, kwargs)

  u0 = get_concrete_u0(prob, tspan[1], kwargs)

  if prob.constant_lags isa Function
    constant_lags = prob.constant_lags(prob.p)
  else
    constant_lags = prob.constant_lags
  end

  u0 = promote_u0(u0, prob.p, tspan[1])
  tspan = promote_tspan(u0, prob.p, tspan, prob, kwargs)

  remake(prob; u0 = u0, tspan = tspan, constant_lags = constant_lags)
end

function get_concrete_tspan(prob, kwargs)
  if prob.tspan isa Function
    tspan = prob.tspan(prob.p)
  elseif prob.tspan === (nothing, nothing)
    if haskey(kwargs, :tspan)
      tspan = kwargs[:tspan]
    else
      error("No tspan is set in the problem or chosen in the init/solve call")
    end
  else
    tspan = prob.tspan
  end

  tspan
end

function get_concrete_u0(prob, t0, kwargs)
  if eval_u0(prob.u0)
    u0 = prob.u0(prob.p, t0)
  elseif prob.u0 === nothing
    u0 = kwargs[:u0]
  else
    u0 = prob.u0
  end

  handle_distribution_u0(u0)
end

handle_distribution_u0(_u0) = _u0
eval_u0(u0::Function) = true
eval_u0(u0) = false

"""
$(SIGNATURES)

Check whether the values of `u0` and `tspan` are appropriate for use with
adaptive integrators and emit specific warnings if they are not.
"""
function adaptive_warn(u0,tspan)
  adaptive_integer_warn(tspan)
end

"""
$(SIGNATURES)

Emit a warning about incompatibility with adaptive integers if `tspan` contains
integers.
"""
function adaptive_integer_warn(tspan)
  if eltype(tspan) <: Integer
    @warn("Integer time values are incompatible with adaptive integrators. Utilize floating point numbers instead of integers in this case, i.e. (0.0,1.0) instead of (0,1).")
  end
end

function __solve(prob::DEProblem,args...;default_set=false,second_time=false,kwargs...)
  if second_time
    error("Default algorithm choices require DifferentialEquations.jl. Please specify an algorithm or import DifferentialEquations directly.")
  elseif length(args) > 0 && !(typeof(args[1]) <: Union{Nothing,DEAlgorithm})
    error("Inappropiate solve command. The arguments do not make sense. Likely, you gave an algorithm which does not actually exist (or does not `<:DiffEqBase.DEAlgorithm`)")
  else
    __solve(prob::DEProblem,nothing,args...;default_set=false,second_time=true,kwargs...)
  end
end

################### Concrete Solve

function _concrete_solve end

function concrete_solve(prob::DiffEqBase.DEProblem,alg::DiffEqBase.DEAlgorithm,
                        u0=prob.u0,p=prob.p,args...;kwargs...)
  _concrete_solve(prob,alg,u0,p,args...;kwargs...)
end

function _concrete_solve(prob::DiffEqBase.DEProblem,alg::DiffEqBase.DEAlgorithm,
                        u0=prob.u0,p=prob.p,args...;kwargs...)
  sol = solve(remake(prob,u0=u0,p=p),alg,args...;kwargs...)
  RecursiveArrayTools.DiffEqArray(sol.u,sol.t)
end

function DiffEqBase._concrete_solve(prob::DiffEqBase.AbstractSteadyStateProblem,
  alg::DiffEqBase.DEAlgorithm, u0 = prob.u0, p = prob.p, args...; kwargs...)
  solve(remake(prob, u0 = u0, p = p), alg, args...; kwargs...)
end

function ChainRulesCore.frule(::typeof(concrete_solve),prob,alg,u0,p,args...;
                     sensealg=nothing,kwargs...)
  _concrete_solve_forward(prob,alg,sensealg,u0,p,args...;kwargs...)
end

function ChainRulesCore.rrule(::typeof(concrete_solve),prob,alg,u0,p,args...;
                     sensealg=nothing,save_idxs=nothing,kwargs...)
  y, back = _concrete_solve_adjoint(prob,alg,sensealg,u0,p,args...;kwargs...)
  isnothing(save_idxs) && return y, back
  y[save_idxs, :], Δ -> (Δ′ = zero(y); Δ′[save_idxs, :] .= Δ; back(Δ′))
end

ZygoteRules.@adjoint function concrete_solve(prob,alg,u0,p,args...;
                      sensealg=nothing,save_idxs=nothing,kwargs...)
  y, back = _concrete_solve_adjoint(prob,alg,sensealg,u0,p,args...;kwargs...)
  isnothing(save_idxs) && return y, back
  y[save_idxs, :], Δ -> (Δ′ = zero(y); Δ′[save_idxs, :] .= Δ; back(Δ′))
end

function _concrete_solve_adjoint(args...;kwargs...)
  error("No adjoint rules exist. Check that you added `using DiffEqSensitivity`")
end

function _concrete_solve_forward(args...;kwargs...)
  error("No sensitivity rules exist. Check that you added `using DiffEqSensitivity`")
end