text:
### The following code defines all the policies that can be used in the simulator

### A large family of policies are priority policies (based on a priority rule
### such as order of arrival, deadline order, etc)
### The general_priority_policty implements this general case. perm is the priority permutation

function general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})

    p=0.0;
    i=1;

    while p<C && i<=length(evs)
        alloc = min(evs[perm[i]].chargingPower,C-p);
        p=p+alloc;
        evs[perm[i]].currentPower=alloc;
        i=i+1;
    end

    while i<=length(evs)
        evs[perm[i]].currentPower=0.0;
        i=i+1
    end

    return p;
end

### Priority based policies

### Earliest Deadline first.
function edf_policy(evs::Array{EVinstance},C::Number)

    deadlines = [ev.currentReportedDeadline for ev in evs];
    perm = sortperm(deadlines);

    return general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})
end

@addpolicy("edf")

### Least laxity first
function llf_policy(evs::Array{EVinstance},C::Number)

    laxities = [ev.currentReportedDeadline-ev.currentWorkload/ev.chargingPower for ev in evs];
    perm = sortperm(laxities);

    return general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})
end

@addpolicy("llf")

function llr_policy(evs::Array{EVinstance},C::Number)

    relative_laxities = [ev.currentReportedDeadline*ev.chargingPower/ev.currentWorkload for ev in evs];
    perm = sortperm(relative_laxities);

    return general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})
end

@addpolicy("llr")

### FIFO. It's a priority policy with arrival time as priority vector.
function fifo_policy(evs::Array{EVinstance},C::Number)

    perm = 1:length(evs)

    return general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})
end

@addpolicy("fifo")

### FIFO. It's a priority policy with reversed arrival time as priority vector.
function lifo_policy(evs::Array{EVinstance},C::Number)

    perm = length(evs):-1:1

    return general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})
end

@addpolicy("lifo")

### LAR: least attained ratio
function lar_policy(evs::Array{EVinstance},C::Number)

    relative_attained = [(ev.departureTime-ev.arrivalTime-ev.currentReportedDeadline)*ev.chargingPower/ev.currentWorkload for ev in evs];
    perm = sortperm(relative_attained,rev=true);

    return general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})
end

@addpolicy("lar")

### LAS: least attained service
function las_policy(evs::Array{EVinstance},C::Number)

    attained = [ev.requestedEnergy-ev.currentWorkload for ev in evs];
    perm = sortperm(attained);

    return general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})
end

@addpolicy("las")

### RATIO policy: just sort by current percentage of charge remaining
function ratio_policy(evs::Array{EVinstance},C::Number)

    ratios = [ev.currentWorkload/ev.requestedEnergy for ev in evs];
    perm = sortperm(ratios,rev=true);

    return general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})
end

@addpolicy("ratio")

### LRPT: Largest remaining processing tiem
function lrpt_policy(evs::Array{EVinstance},C::Number)

    remaining = [ev.currentWorkload for ev in evs];
    perm = sortperm(remaining,rev=true);

    return general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})
end

@addpolicy("lrpt")

### Non-priority based policies

### Processor sharing Policy
function parallel_policy(evs::Array{EVinstance},C::Number)

    totPower=sum([ev.chargingPower for ev in evs]);
    curtail = min(1,C/totPower);

    [ev.currentPower = curtail*ev.chargingPower for ev in evs]

    return curtail*totPower
end

@addpolicy("parallel")

### Proportional fairness policy
function pf_policy(evs::Array{EVinstance},C::Number)

    ##TODO Revisar comput_pf. Parece estar bien.
    workloads = [ev.currentWorkload for ev in evs];
    deadlines = [ev.currentReportedDeadline for ev in evs];
    U=compute_pf(workloads,deadlines,C)

    for i=1:length(evs)
        evs[i].currentPower = U[i]
    end

    return sum(U);
end

function compute_pf(workloads,deadlines,C)

    w=workloads./deadlines;

    U=zeros(size(w));
    perm = sortperm(w,rev=true);
    w=sort(w,rev=true);

    for i=1:length(w)

        aux = w[i:end]/sum(w[i:end])*C;
        if aux[1]>1
            U[i] = 1;
            C=C-1;
        else
            U[i:end] = aux;
            break
        end
    end
    U[perm]=U;
    return U;

end

@addpolicy("pf")

### Exact scheduling.
function exact_policy(evs::Array{EVinstance},C::Number)

    #exact scheduling o potencia maxima
    U = [min(ev.currentWorkload/ev.currentReportedDeadline,ev.chargingPower) for ev in evs];

    #curtailing si me paso de C
    if sum(U)>C
        U = C/sum(U)*U;
        #otra posibilidad es ordernar por rate y asignar hasta C
    end

    for i=1:length(evs)
        evs[i].currentPower = U[i]
    end

    return sum(U);
end

@addpolicy("exact")

### MPC peak minimizer policy. Check whether this is needed to reduce dependencies.
function peak_policy(evs::Array{EVinstance},C::Number)

    U=zeros(length(evs));

    idx = sortperm([ev.currentReportedDeadline for ev in evs]);

    sigma = [ev.currentWorkload for ev in evs][idx];
    tau = [ev.currentReportedDeadline for ev in evs][idx];
    deltat=diff([0;tau]);
    power = [ev.chargingPower for ev in evs][idx];
    n=length(evs);

    m=Model(GLPK.Optimizer)

    @variable(m,x[1:n,1:n]>=0)
    @variable(m,auxvar)

    @constraint(m,[i=1:n,j=i+1:n],x[i,j]==0)

    @constraint(m,[i=1:n,j=1:i],x[i,j]<=power[i])
    @constraint(m,[i=1:n],sum(x[i,:].*deltat)==sigma[i])
    @constraint(m,sum(x,dims=1).<=auxvar)

    @objective(m,Min,auxvar)

    solve(m)

    U[idx] = max.(getvalue(x)[:,1],0.0);

    for i=1:length(evs)
        evs[i].currentPower = U[i]
    end

    return sum(U);

end

@addpolicy("peak")

#### WEIRD POLICIES: This are only for trial purposes. Will be depurated.

#max weight policy where weight is minimum between rem. work and rem. deadline
function mw_policy(evs::Array{EVinstance},C::Number)

    remaining_w = [ev.currentWorkload for ev in evs];
    remaining_d = [ev.currentReportedDeadline for ev in evs];

    weights = min.(remaining_w,remaining_d)
    perm = sortperm(weights,rev=true);

    return general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})
end

@addpolicy("mw")

#this policy comes from maximizing myopically the potential amount of work one can perform, ignoring future arrivals
#it underload it behaves exactly as exact scheduling!
function weird_policy(evs::Array{EVinstance},C::Number)

    perm = sortperm([ev.currentReportedDeadline for ev in evs], rev=true)

    remaining_w = [ev.currentWorkload for ev in evs];
    remaining_d = [ev.currentReportedDeadline for ev in evs];

    return general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})
end

@addpolicy("weird")

### Curtailed policies (for reported deadlines)
function edfc_policy(evs::Array{EVinstance},C::Number)

    deadlines = [ev.currentReportedDeadline for ev in evs];

    positivas = findall(deadlines.>0)
    negativas = findall(deadlines.<=0)

    perm1 = sortperm(deadlines[positivas]);
    perm2 = sortperm(deadlines[negativas], rev=true);

    perm = [positivas[perm1];negativas[perm2]]

    return general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})

end

@addpolicy("edfc")

function llfc_policy(evs::Array{EVinstance},C::Number)

    deadlines = [ev.currentReportedDeadline for ev in evs];
    laxities = [ev.currentReportedDeadline-ev.currentWorkload/ev.chargingPower for ev in evs];

    positivas = findall(deadlines.>0)
    negativas = findall(deadlines.<=0)

    perm1 = sortperm(laxities[positivas]);
    perm2 = sortperm(laxities[negativas], rev=true);

    perm = [positivas[perm1];negativas[perm2]]

    return general_priority_policy(evs::Array{EVinstance},C::Number,perm::Array{<:Integer})

end

@addpolicy("llfc")


text:
using ActuaryUtilities

using Dates
using Test

const Yields = ActuaryUtilities.Yields

import DayCounts

include("risk_measures.jl")

@testset "Temporal functions" begin
    @testset "years_between" begin
        @test years_between(Date(2018, 9, 30), Date(2018, 9, 30)) == 0
        @test years_between(Date(2018, 9, 30), Date(2018, 9, 30), true) == 0
        @test years_between(Date(2018, 9, 30), Date(2019, 9, 30), false) == 0
        @test years_between(Date(2018, 9, 30), Date(2019, 9, 30), true) == 1
        @test years_between(Date(2018, 9, 30), Date(2019, 10, 1), true) == 1
        @test years_between(Date(2018, 9, 30), Date(2019, 10, 1), false) == 1
    end

    @testset "duration tests" begin
        @test duration(Date(2018, 9, 30), Date(2019, 9, 30)) == 2
        @test duration(Date(2018, 9, 30), Date(2018, 9, 30)) == 1
        @test duration(Date(2018, 9, 30), Date(2018, 10, 1)) == 1
        @test duration(Date(2018, 9, 30), Date(2019, 10, 1)) == 2
        @test duration(Date(2018, 9, 30), Date(2018, 6, 30)) == 0
        @test duration(Date(2018, 9, 30), Date(2017, 6, 30)) == -1
        @test duration(Date(2018, 10, 15), Date(2019, 9, 30)) == 1
        @test duration(Date(2018, 10, 15), Date(2019, 10, 30)) == 2
        @test duration(Date(2018, 10, 15), Date(2019, 10, 15)) == 2
        @test duration(Date(2018, 10, 15), Date(2019, 10, 14)) == 1
    end
end

@testset "accum_offset" begin
    @test all(accum_offset([0.9, 0.8, 0.7]) .== [1.0,0.9,1.0 * 0.9 * 0.8])
    @test all(accum_offset([0.9, 0.8, 0.7],op=+) .== [1.0,1.9,2.7])
    @test all(accum_offset([0.9, 0.8, 0.7],op=+,init=2) .== [2.0,2.9,3.7])

    @test all(accum_offset([1, 2, 3]) .== [1,1,2])
end

@testset "financial calcs" begin

    @testset "pv" begin
        cf = [100, 100]

        @test pv(0.05, cf) ≈ cf[1] / 1.05 + cf[2] / 1.05^2
        @test price(0.05, cf) ≈ pv(0.05, cf)

        # this vector came from Numpy Financial's test suite with target of 122.89, but that assumes payments are begin of period
        # 117.04 comes from Excel verification with NPV function
        @test isapprox(pv(0.05, [-15000, 1500, 2500, 3500, 4500, 6000]), 117.04, atol = 1e-2)

        cfs = ones(3)
        @test present_values(Yields.Constant(0.0),cfs) == [3,2,1]
        pvs = present_values(Yields.Constant(0.1),cfs)
        @test pvs[3] ≈ 1 / 1.1
        @test pvs[2] ≈ (1 / 1.1 + 1) / 1.1

        @test all(present_values(0.00, [1,1,1]) .≈ [3,2,1])
        @test all(present_values(0.00, [1,1,1],[0,1,2]) .≈ [3,2,1])
        @test all(present_values(0.00, [1,1,1],[1,2,3]) .≈ [3,2,1])
        @test all(present_values(0.00, [1,1,1],[1,2,3]) .≈ [3,2,1])
        @test all(present_values(0.01, [1,2,3]) .≈ [ 5.862461552497766,4.921086168022744,2.9702970297029707])
        @test all(present_values(Yields.Forward([0.1,0.2]), [10,20],[0,1]) ≈  [28.18181818181818,18.18181818181818])
        @test all(present_values([0.1,0.2], [10,20],[0,1]) ≈  [28.18181818181818,18.18181818181818])
    end

    @testset "pv with timepoints" begin
        cf = [100, 100]

        @test pv(0.05, cf, [1,2]) ≈ cf[1] / 1.05 + cf[2] / 1.05^2
    end

    @testset "pv with vector discount rates" begin
        cf = [100, 100]
        @test pv([0.0,0.05], cf) ≈ 100 / 1.0 + 100 / 1.05
        @test pv(ActuaryUtilities.Yields.Forward([0.0,0.05]), cf) ≈ 100 / 1.0 + 100 / 1.05
        @test pv([0.05,0.0], cf) ≈ 100 / 1.05 + 100 / 1.05
        @test pv([0.05,0.1], cf) ≈ 100 / 1.05 + 100 / 1.05 / 1.1

        ts = [0.5,1]
        @test pv(ActuaryUtilities.Yields.Forward([0.0,0.05], ts), cf, ts) ≈ 100 / 1.0 + 100 / 1.05^0.5
        @test pv(ActuaryUtilities.Yields.Forward([0.05,0.0], ts), cf, ts) ≈ 100 / 1.05^0.5 + 100 / 1.05^0.5
        @test pv(ActuaryUtilities.Yields.Forward([0.05,0.1], ts), cf, ts) ≈ 100 / 1.05^0.5 + 100 / (1.05^0.5) / (1.1^0.5)

        #without explicit Yields constructor
        @test pv([0.0,0.05], cf, ts) ≈ 100 / 1.0 + 100 / 1.05^0.5

        @test price([0.0,0.05], cf, ts) ≈ pv([0.0,0.05], cf, ts)
        @test price([0.0,0.05], -1 .* cf, ts) ≈ abs(pv([0.0,0.05], cf, ts))

    end

    @testset "irr" begin

        v = [-70000,12000,15000,18000,21000,26000]

        # per Excel (example comes from Excel help text)
        @test isapprox(irr(v[1:2]), -0.8285714285714, atol = 0.001)
        @test isapprox(ActuaryUtilities.irr(v[1:2]), -0.8285714285714, atol = 0.001)
        @test isapprox(ActuaryUtilities.irr(v[1:3]), -0.4435069413346, atol = 0.001)
        @test isapprox(irr(v[1:4]), -0.1821374641455, atol = 0.001)
        @test isapprox(irr(v[1:5]), -0.0212448482734, atol = 0.001)
        @test isapprox(irr(v[1:6]),  0.0866309480365, atol = 0.001)

        @test_throws MethodError irr("hello")

        # much more challenging to solve b/c of the overflow below zero
        cfs = [t % 10 == 0 ? -10 : 1.5 for t in 0:99]

        @test isapprox(irr(cfs), 0.06463163963925866, atol = 0.001)

        # issue #28
        cfs = [-8.728037307132952e7, 3.043754023830998e7, 2.963004184784189e7, 2.8803030748755097e7, 2.7956912111811966e7, 2.7092182051244527e7, 2.6209069543806538e7, 2.5307964329840004e7, 2.438961041057478e7, 2.3455084653011695e7, 2.2505925520018265e7, 2.154395414765592e7, 2.0571076113065004e7, 1.958930608135183e7, 1.8600627464895025e7, 1.7606980923262402e7, 1.661046149512893e7, 1.561312825963898e7, 1.461760481586352e7, 1.3626801207410209e7, 1.2644733969499402e7, 1.1675393687299855e7, 1.0722720151658386e7, 9.79075673433771e6, 8.883278741880089e6, 8.004445298876338e6, 7.1588010859461725e6, 6.351121678665243e6, 5.585860320479795e6, 4.8673895159943625e6, 4.19908059495347e6, 3.583538247530099e6, 3.022766488834396e6, 2.5181072324190177e6, 2.0701053881076649e6, 1.6782921224664208e6, 1.3410605489291362e6, 1.0556643097527474e6, 818348.5357315112, 624147.9373214925, 467849.788997191, 344241.752520618, 248285.65630649775, 175235.5475426321, 120677.87174498942, 80759.09804678289, 52186.83400936739, 32211.057718402008, 18589.51907385164, 9540.782278174447, 3688.4015341755294]
        @test irr(cfs,0:50) ≈ 0.3176680627111823

        @test irr([-100,100]) ≈ 0.
        @test isnothing(irr([100,100])) # answer is -1, but search range won't find it

        # test the unsolvable
        @test isnothing(irr([-1e8,0.,0.,0.],0:3))

    end

    @testset "irr with fractional time" begin
        irr1 = irr([-10,5,5,5],[0,1,2,3])
        @test irr1 ≈ irr([-10,5,5,5])
        irr2 = irr([-10,5,5,5],[0,1,2,3] ./ 2)

        @test (1+irr1)^2-1 ≈ irr2

    end

    @testset "numpy examples" begin

        @test isapprox(irr([-150000, 15000, 25000, 35000, 45000, 60000]),  0.0524,     atol = 1e-4)
        @test isapprox(irr([-100, 0, 0, 74]), -0.0955,     atol = 1e-4)
        @test isapprox(irr([-100, 39, 59, 55, 20]),  0.28095,    atol = 1e-4)
        @test isapprox(irr([-100, 100, 0, -7]), -0.0833,     atol = 1e-4)
        @test isapprox(irr([-100, 100, 0, 7]),  0.06206,    atol = 1e-4)

        # this has multiple roots, of which 0.709559 and 0.0886. Want to find the one closer to zero
        @test isapprox(irr([-5, 10.5, 1, -8, 1]),  0.0886,     atol = 1e-4)
    end

    @testset "xirr with float times" begin

        @test isapprox(irr([-100,100], [0,1]), 0.0, atol = 0.001)
        @test isapprox(irr([-100,110], [0,1]), 0.1, atol = 0.001)

    end

    @testset "xirr with real dates" begin

        v = [-70000,12000,15000,18000,21000,26000]
        dates = Date(2019, 12, 31):Year(1):Date(2024, 12, 31)
        times = map(d->DayCounts.yearfrac(dates[1], d, DayCounts.Thirty360()), dates)
    # per Excel (example comes from Excel help text)
        @test isapprox(irr(v[1:2], times[1:2]), -0.8285714285714, atol = 0.001)
        @test isapprox(ActuaryUtilities.irr_robust(v[1:3], times[1:3]), -0.4435069413346, atol = 0.001)
        @test isapprox(irr(v[1:4], times[1:4]), -0.1821374641455, atol = 0.001)
        @test isapprox(irr(v[1:5], times[1:5]), -0.0212448482734, atol = 0.001)
        @test isapprox(irr(v[1:6], times[1:6]),  0.0866309480365, atol = 0.001)

    end
end

@testset "Breakeven time" begin

    @testset "basic" begin
        @test breakeven(0.10, [-10,1,2,3,4,8]) == 5
        @test breakeven(0.10, [-10,15,2,3,4,8]) == 1
        @test breakeven(0.10, [-10,15,2,3,4,8]) == 1
        @test breakeven(0.10, [10,15,2,3,4,8]) == 0
        @test isnothing(breakeven(0.10, [-10,-15,2,3,4,8]))
    end

    @testset "basic with vector interest" begin
        @test breakeven(0.0,[-10,1,2,3,4], [1,2,3,4,5]) == 5
        #
        @test isnothing(breakeven([0.0,0.0,0.0,0.0,0.1], [-10,1,2,3,4], [1,2,3,4,5]))
        @test breakeven([0.0,0.0,0.0,0.0,-0.5], [-10,1,2,3,4], [1,2,3,4,5]) == 5
        @test breakeven([0.0,0.0,0.0,-0.9,-0.5],[-10,1,2,3,4], [1,2,3,4,5]) == 4
        @test breakeven([0.1,0.1,0.2,0.1,0.1], [-10,1,12,3,4], [1,2,3,4,5]) == 3
    end

    @testset "timepoints" begin
        times = [t for t in 0:5]
        @test breakeven(0.10,[-10,1,2,3,4,8], times) == 5
        @test breakeven(0.10,[-10,15,2,3,4,8], times) == 1
        @test breakeven(0.10,[-10,15,2,3,4,8], times) == 1
        @test isnothing(breakeven(0.10,[-10,-15,2,3,4,8], times))
    end
end

@testset "duration and convexity" begin

    @testset "wikipedia example" begin
        times = [0.5,1,1.5,2]
        cfs = [10,10,10,110]
        V = present_value(0.04, cfs, times)

        @test duration(Macaulay(), 0.04, cfs, times) ≈ 1.777570320376649
        @test duration(Modified(), 0.04, cfs, times) ≈ 1.777570320376649 / (1 + 0.04)
        @test duration(0.04, cfs, times) ≈ 1.777570320376649 / (1 + 0.04)
        @test duration(DV01(), 0.04, cfs, times) ≈ 1.777570320376649 / (1 + 0.04) * V / 100

        @test duration(Yields.Constant(0.04), cfs, times) ≈ 1.777570320376649 / (1 + 0.04)
        @test duration(Yields.Constant(0.04), -1 .* cfs, times) ≈ 1.777570320376649 / (1 + 0.04) atol=0.00001
        @test duration(Yields.Forward([0.04,0.04]), cfs, times) ≈ 1.777570320376649 / (1 + 0.04) atol=0.00001

        # test that dispatch resolves the ambiguity between duration(Yield,vec) and duration(Yield, function)
        @test duration(Yields.Constant(0.03),cfs) > 0
        @test convexity(Yields.Constant(0.03),cfs) > 0
    end

    @testset "finpipe example" begin
        # from https://www.finpipe.com/duration-macaulay-and-modified-duration-convexity/

        cfs = zeros(10) .+ 3.75
        cfs[10] += 100

        times = 0.5:0.5:5.0
        int = (1 + 0.075 / 2)^2 - 1 # convert bond yield to effective yield

        @test isapprox(present_value(int, cfs, times), 100.00, atol = 1e-2)
        @test isapprox(duration(Macaulay(), int, cfs, times), 4.26, atol = 1e-2)
    end

    @testset "Primer example" begin
        # from https://math.illinoisstate.edu/krzysio/Primer.pdf
        # the duration tests are commented out because I think the paper is wrong on the duration?
        cfs = [0,0,0,0,1.0e6]
        times = 1:5

        @test isapprox(present_value(0.04, cfs, times), 821927.11, atol = 1e-2)
        # @test isapprox(duration(0.04,cfs,times),4.76190476,atol=1e-6)
        @test isapprox(convexity(0.04, cfs, times), 27.7366864, atol = 1e-6)
        @test isapprox(convexity(0.04, cfs), 27.7366864, atol = 1e-6)
        # the same, but with a functional argument
        value(i) = present_value(i, cfs, times)
        # @test isapprox(duration(0.04,value),4.76190476,atol=1e-6)
        @test isapprox(convexity(0.04, value), 27.7366864, atol = 1e-6)
    end

    @testset "Quantlib" begin
    # https://mhittesdorf.wordpress.com/2013/03/12/introduction-to-quantlib-duration-and-convexity/
        cfs = [5,5,105]
        times = 1:3
        @test present_value(0.03, cfs, times) ≈ 105.6572227097894
        @test duration(Macaulay(), 0.03, cfs, times) ≈ 2.863504670671131
        @test duration(0.03, cfs, times) ≈ 2.780101622010806
        @test convexity(0.03, cfs, times) ≈ 10.62580548268594

        # test omitting the times argument
        @test duration(Macaulay(), 0.03, cfs) ≈ 2.863504670671131
        @test duration(0.03, cfs) ≈ 2.780101622010806
        @test convexity(0.03, cfs) ≈ 10.62580548268594

        # test a single matrix dimension
        cfs = [5 0 0
               0 5 105]

        @test duration(0.03, sum(cfs,dims=1), times) ≈ 2.780101622010806

        cfs = [5 0
               5 0
               0 105]

        @test duration(0.03, sum(cfs,dims=2), times) ≈ 2.780101622010806

    end

end

include("run_doctests.jl")

ext: jl


text:
function t_test(n,r)
    stat = r *sqrt((n-2)/1-r^2)
    t_dist = TDist(n-2)
    return pvalue(t_dist,stat,tail=:right)
end

function rank_correlation_test(radiance, clouds)
    missings = findall(ismissing, radiance)
    y = filter!(!ismissing, copy(radiance))
    x = Array{Union{Float64, Missing}}(copy(clouds))
    x[missings] .= missing
    x = Array{Float64}(filter!(!ismissing, x))
    y_de = Array{Float64}(detrend_ts(y))
    r = corspearman(y_de, x)
    n = length(y_de)
    return t_test(n,r)
end


text:
using Combinatorics
lines = readlines("inputs/6-10/08.txt")

function part1(lines)
  t = 0
  for line in lines
    lin = split(line, ' ')[end-4:end]
    foo = filter(x -> length(x) ∈ [2, 3, 4, 7], lin)
    t += length(foo)
  end
  t
end

str_sort(str) = str |> collect |> sort |> join

# brute force solution that I submitted
function part2slow(lines)
  d = Dict("cf" => 1,
    "acdeg" => 2,
    "acdfg" => 3,
    "bcdf" => 4,
    "abdfg" => 5,
    "abdefg" => 6,
    "acf" => 7,
    "abcdefg" => 8,
    "abcdfg" => 9,
    "abcefg" => 0)

  score(new_ans, d) = parse(Int, join([d[join(str_sort(x))] for x in new_ans]))

  t = 0
  for line in lines
    wires, ans = split(line, '|')
    wires = split(wires, " ")[1:end-1]
    ans = split(ans, " ")[2:end]
    for perm in permutations("abcdefg" |> collect)
      mapping = Dict(i => j for (i, j) in zip(perm, "abcdefg"))
      new_wires = map(x -> join(mapping[f] for f in x), wires)
      new_ans = map(x -> join(mapping[f] for f in x), ans)

      if all((join(str_sort(wire))) in keys(d) for wire in new_wires)
        t += score(new_ans, d)
        break
      end
    end
  end
  t
end

score(new_ans, d) = parse(Int, join([d[join(str_sort(x))] for x in new_ans]))

lines = readlines("inputs/6-10/08.txt")

strdiff(a, b) = setdiff(Set(a), Set(b)) |> collect

num_in(char, str) = count(x -> x == char, collect(str))

d = Dict("cf" => 1,
  "acdeg" => 2,
  "acdfg" => 3,
  "bcdf" => 4,
  "abdfg" => 5,
  "abdefg" => 6,
  "acf" => 7,
  "abcdefg" => 8,
  "abcdfg" => 9,
  "abcefg" => 0)

function solve_line(line)
  alphabet = "abcdefg"
  s = split(line, ' ')
  wires = s[1:end-5]
  ans = s[end-3:end]
  nums = Dict()
  segs = Dict()
  charcounts = Dict()
  nums[1] = wires[findfirst(x -> length(x) == 2, wires)]
  nums[7] = wires[findfirst(x -> length(x) == 3, wires)]
  nums[4] = wires[findfirst(x -> length(x) == 4, wires)]
  nums[8] = wires[findfirst(x -> length(x) == 7, wires)]
  segs['a'] = strdiff(nums[7], nums[1])[1]
  activated_wires = join(wires) |> str_sort
  for char in collect(alphabet)
    charcounts[char] = num_in(char, activated_wires)
    activations = num_in(char, activated_wires)
    if activations == 6
      segs['b'] = char
    elseif activations == 9
      segs['f'] = char
    elseif activations == 4
      segs['e'] = char
    elseif activations == 8 && char != segs['a']
      segs['c'] = char
    end
  end
  segs['d'] = strdiff(nums[4], join(values(segs)))[1]
  segs['g'] = strdiff(alphabet, join(values(segs)))[1]

  rev = Dict(value => key for (key, value) in segs)
  parse(Int, [d[x] for x in map(an -> [rev[x] for x in an] |> join |> str_sort, ans)] |> join)
end

function part2(lines)
  t = 0
  for line in lines
    t += solve_line(line)
  end
  t
end

@show part1(lines)
@show part2(lines)


text:
tV = SimpleRDArray(V)
tM = SimpleRDArray(M)
tT = SimpleRDArray(T)

@testset "Unary Operations: $f" for f in (:+, :-, :conj, :real, :imag)
    @eval begin
        @test $f(tV) == $f(V)
        @test $f(tM) == $f(M)
        @test $f(tT) == $f(T)
    end
end

@testset "Binary Operations: $f" for f in (:+, :-)
    @eval begin
        @test $f(tV, tV) == $f(tV, V) == $f(V, tV) == $f(V, V)
        @test $f(tM, tM) == $f(tM, M) == $f(M, tM) == $f(M, M)
        @test $f(tT, tT) == $f(tT, T) == $f(T, tT) == $f(T, T)
    end
end

@testset "N-ary Operations: $f" for f in (:+,)
    @eval begin
        @test $f(tV, tV, tV) == $f(V, V, V)
        @test $f(tM, tM, tM) == $f(M, M, M)
        @test $f(tT, tT, tT) == $f(T, T, T)
    end
end

@testset "Broadcast: $f" for f in (:+, :-, :*, :/, :\, :^)
    @eval begin
        @test $f.(tV, tV) == $f.(tV, V) == $f.(V, tV) == $f.(V, V)
    end
end

@testset "Linear Algebra: $f" for f in (transpose, adjoint)
    @eval @test $f(tV) == $f(V)
    @eval @test $f(tM) == $f(M)
    for bf in (:+, :-)
        @eval @test $bf($f(tV), $f(tV)) == $bf($f(V), $f(V))
        @eval @test $bf($f(tM), $f(tM)) == $bf($f(M), $f(M))
    end
    for bf in (:*,)
        @eval begin
            @test $bf($f(tV), tV) == $bf($f(V), V)
            @test $bf(tV, $f(tV)) == $bf(V, $f(V))
            @test $bf($f(tM), tM) == $bf($f(M), M)
            @test $bf(tM, $f(tM)) == $bf(M, $f(M))
        end
    end
end

@testset "Linear Algebra: M $f V and M $f M" for f in (:*, :\)
    @eval begin
        @test $f(tM, tV) == $f(tM, V) == $f(M, tV) == $f(M, V)
        @test $f(tM, tM) == $f(tM, M) == $f(M, tM) == $f(M, M)
    end
end


text:
module Day13

using Advent: save_input
using Mods: CRT, Mod, modulus

function solve_a(filename)
    open(filename) do io
        timestamp = parse(Int, readline(io))
        times = [parse(Int, x) for x in split(readline(io), ",") if x != "x"]

        minutes, i = findmin(times .- timestamp .% times)
        minutes * times[i]
    end
end

function solve_b(filename)
    line = readlines(filename)[2]
    times = (tryparse(Int, x) for x in split(line, ","))
    mods = (Mod(s - 1, t) for (s, t) in enumerate(times) if t !== nothing)
    res = CRT(mods...)
    modulus(res) - res.val
end

solve(filename = save_input(2020, 13)) = solve_a(filename), solve_b(filename)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using DelimitedFiles

# Converting V_0 and wfc into formats easily plotted by gnuplot
function convert_to_gnuplot(xDim, wfc_file, wfci_file, outfile)
    wfc = reshape(readdlm(wfc_file) + readdlm(wfci_file)*im,256, 256)

    writedlm(outfile, abs2.(wfc[:,div(xDim,2)]))
end

data_dir = "./"
convert_to_gnuplot(256, data_dir*"wfc_g0_0",data_dir*"wfc_g0i_0",
                   data_dir*"wfc_g0.dat")
convert_to_gnuplot(256, data_dir*"wfc_g0.5_0",data_dir*"wfc_g0.5i_0",
                   data_dir*"wfc_g0.5.dat")
convert_to_gnuplot(256, data_dir*"wfc_g1_0",data_dir*"wfc_g1i_0",
                   data_dir*"wfc_g1.dat")
convert_to_gnuplot(256, data_dir*"wfc_g1.5_0",data_dir*"wfc_g1.5i_0",
                   data_dir*"wfc_g1.5.dat")
convert_to_gnuplot(256, data_dir*"wfc_g2_0",data_dir*"wfc_g2i_0",
                   data_dir*"wfc_g2.dat")

max_stars_repo_name: leios/thesis

max_issues_repo_name: leios/thesis

max_forks_repo_name: leios/thesis


text:
# Julia wrapper for header: /Users/jgoldfar/Documents/misc/julia/Hypre/deps/usr/include/HYPRE_lobpcg.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function HYPRE_AssumedPartitionCheck()
    ccall((:HYPRE_AssumedPartitionCheck,libHypre),HYPRE_Int,())
end

function HYPRE_PCGSetup(solver::HYPRE_Solver,A::HYPRE_Matrix,b::HYPRE_Vector,x::HYPRE_Vector)
    ccall((:HYPRE_PCGSetup,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector,HYPRE_Vector),solver,A,b,x)
end

function HYPRE_PCGSolve(solver::HYPRE_Solver,A::HYPRE_Matrix,b::HYPRE_Vector,x::HYPRE_Vector)
    ccall((:HYPRE_PCGSolve,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector,HYPRE_Vector),solver,A,b,x)
end

function HYPRE_PCGSetTol(solver::HYPRE_Solver,tol::HYPRE_Real)
    ccall((:HYPRE_PCGSetTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,tol)
end

function HYPRE_PCGSetAbsoluteTol(solver::HYPRE_Solver,a_tol::HYPRE_Real)
    ccall((:HYPRE_PCGSetAbsoluteTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,a_tol)
end

function HYPRE_PCGSetResidualTol(solver::HYPRE_Solver,rtol::HYPRE_Real)
    ccall((:HYPRE_PCGSetResidualTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,rtol)
end

function HYPRE_PCGSetAbsoluteTolFactor(solver::HYPRE_Solver,abstolf::HYPRE_Real)
    ccall((:HYPRE_PCGSetAbsoluteTolFactor,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,abstolf)
end

function HYPRE_PCGSetConvergenceFactorTol(solver::HYPRE_Solver,cf_tol::HYPRE_Real)
    ccall((:HYPRE_PCGSetConvergenceFactorTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,cf_tol)
end

function HYPRE_PCGSetStopCrit(solver::HYPRE_Solver,stop_crit::HYPRE_Int)
    ccall((:HYPRE_PCGSetStopCrit,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,stop_crit)
end

function HYPRE_PCGSetMaxIter(solver::HYPRE_Solver,max_iter::HYPRE_Int)
    ccall((:HYPRE_PCGSetMaxIter,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,max_iter)
end

function HYPRE_PCGSetTwoNorm(solver::HYPRE_Solver,two_norm::HYPRE_Int)
    ccall((:HYPRE_PCGSetTwoNorm,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,two_norm)
end

function HYPRE_PCGSetRelChange(solver::HYPRE_Solver,rel_change::HYPRE_Int)
    ccall((:HYPRE_PCGSetRelChange,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,rel_change)
end

function HYPRE_PCGSetRecomputeResidual(solver::HYPRE_Solver,recompute_residual::HYPRE_Int)
    ccall((:HYPRE_PCGSetRecomputeResidual,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,recompute_residual)
end

function HYPRE_PCGSetRecomputeResidualP(solver::HYPRE_Solver,recompute_residual_p::HYPRE_Int)
    ccall((:HYPRE_PCGSetRecomputeResidualP,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,recompute_residual_p)
end

function HYPRE_PCGSetPrecond(solver::HYPRE_Solver,precond::HYPRE_PtrToSolverFcn,precond_setup::HYPRE_PtrToSolverFcn,precond_solver::HYPRE_Solver)
    ccall((:HYPRE_PCGSetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_PtrToSolverFcn,HYPRE_PtrToSolverFcn,HYPRE_Solver),solver,precond,precond_setup,precond_solver)
end

function HYPRE_PCGSetLogging(solver::HYPRE_Solver,logging::HYPRE_Int)
    ccall((:HYPRE_PCGSetLogging,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,logging)
end

function HYPRE_PCGSetPrintLevel(solver::HYPRE_Solver,level::HYPRE_Int)
    ccall((:HYPRE_PCGSetPrintLevel,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,level)
end

function HYPRE_PCGGetNumIterations(solver::HYPRE_Solver,num_iterations)
    ccall((:HYPRE_PCGGetNumIterations,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,num_iterations)
end

function HYPRE_PCGGetFinalRelativeResidualNorm(solver::HYPRE_Solver,norm)
    ccall((:HYPRE_PCGGetFinalRelativeResidualNorm,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,norm)
end

function HYPRE_PCGGetResidual(solver::HYPRE_Solver,residual)
    ccall((:HYPRE_PCGGetResidual,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{Ptr{Void}}),solver,residual)
end

function HYPRE_PCGGetTol(solver::HYPRE_Solver,tol)
    ccall((:HYPRE_PCGGetTol,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,tol)
end

function HYPRE_PCGGetResidualTol(solver::HYPRE_Solver,rtol)
    ccall((:HYPRE_PCGGetResidualTol,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,rtol)
end

function HYPRE_PCGGetAbsoluteTolFactor(solver::HYPRE_Solver,abstolf)
    ccall((:HYPRE_PCGGetAbsoluteTolFactor,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,abstolf)
end

function HYPRE_PCGGetConvergenceFactorTol(solver::HYPRE_Solver,cf_tol)
    ccall((:HYPRE_PCGGetConvergenceFactorTol,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,cf_tol)
end

function HYPRE_PCGGetStopCrit(solver::HYPRE_Solver,stop_crit)
    ccall((:HYPRE_PCGGetStopCrit,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,stop_crit)
end

function HYPRE_PCGGetMaxIter(solver::HYPRE_Solver,max_iter)
    ccall((:HYPRE_PCGGetMaxIter,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,max_iter)
end

function HYPRE_PCGGetTwoNorm(solver::HYPRE_Solver,two_norm)
    ccall((:HYPRE_PCGGetTwoNorm,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,two_norm)
end

function HYPRE_PCGGetRelChange(solver::HYPRE_Solver,rel_change)
    ccall((:HYPRE_PCGGetRelChange,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,rel_change)
end

function HYPRE_GMRESGetSkipRealResidualCheck(solver::HYPRE_Solver,skip_real_r_check)
    ccall((:HYPRE_GMRESGetSkipRealResidualCheck,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,skip_real_r_check)
end

function HYPRE_PCGGetPrecond(solver::HYPRE_Solver,precond_data_ptr)
    ccall((:HYPRE_PCGGetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Solver}),solver,precond_data_ptr)
end

function HYPRE_PCGGetLogging(solver::HYPRE_Solver,level)
    ccall((:HYPRE_PCGGetLogging,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,level)
end

function HYPRE_PCGGetPrintLevel(solver::HYPRE_Solver,level)
    ccall((:HYPRE_PCGGetPrintLevel,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,level)
end

function HYPRE_PCGGetConverged(solver::HYPRE_Solver,converged)
    ccall((:HYPRE_PCGGetConverged,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,converged)
end

function HYPRE_GMRESSetup(solver::HYPRE_Solver,A::HYPRE_Matrix,b::HYPRE_Vector,x::HYPRE_Vector)
    ccall((:HYPRE_GMRESSetup,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector,HYPRE_Vector),solver,A,b,x)
end

function HYPRE_GMRESSolve(solver::HYPRE_Solver,A::HYPRE_Matrix,b::HYPRE_Vector,x::HYPRE_Vector)
    ccall((:HYPRE_GMRESSolve,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector,HYPRE_Vector),solver,A,b,x)
end

function HYPRE_GMRESSetTol(solver::HYPRE_Solver,tol::HYPRE_Real)
    ccall((:HYPRE_GMRESSetTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,tol)
end

function HYPRE_GMRESSetAbsoluteTol(solver::HYPRE_Solver,a_tol::HYPRE_Real)
    ccall((:HYPRE_GMRESSetAbsoluteTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,a_tol)
end

function HYPRE_GMRESSetConvergenceFactorTol(solver::HYPRE_Solver,cf_tol::HYPRE_Real)
    ccall((:HYPRE_GMRESSetConvergenceFactorTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,cf_tol)
end

function HYPRE_GMRESSetStopCrit(solver::HYPRE_Solver,stop_crit::HYPRE_Int)
    ccall((:HYPRE_GMRESSetStopCrit,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,stop_crit)
end

function HYPRE_GMRESSetMinIter(solver::HYPRE_Solver,min_iter::HYPRE_Int)
    ccall((:HYPRE_GMRESSetMinIter,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,min_iter)
end

function HYPRE_GMRESSetMaxIter(solver::HYPRE_Solver,max_iter::HYPRE_Int)
    ccall((:HYPRE_GMRESSetMaxIter,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,max_iter)
end

function HYPRE_GMRESSetKDim(solver::HYPRE_Solver,k_dim::HYPRE_Int)
    ccall((:HYPRE_GMRESSetKDim,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,k_dim)
end

function HYPRE_GMRESSetRelChange(solver::HYPRE_Solver,rel_change::HYPRE_Int)
    ccall((:HYPRE_GMRESSetRelChange,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,rel_change)
end

function HYPRE_GMRESSetSkipRealResidualCheck(solver::HYPRE_Solver,skip_real_r_check::HYPRE_Int)
    ccall((:HYPRE_GMRESSetSkipRealResidualCheck,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,skip_real_r_check)
end

function HYPRE_GMRESSetPrecond(solver::HYPRE_Solver,precond::HYPRE_PtrToSolverFcn,precond_setup::HYPRE_PtrToSolverFcn,precond_solver::HYPRE_Solver)
    ccall((:HYPRE_GMRESSetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_PtrToSolverFcn,HYPRE_PtrToSolverFcn,HYPRE_Solver),solver,precond,precond_setup,precond_solver)
end

function HYPRE_GMRESSetLogging(solver::HYPRE_Solver,logging::HYPRE_Int)
    ccall((:HYPRE_GMRESSetLogging,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,logging)
end

function HYPRE_GMRESSetPrintLevel(solver::HYPRE_Solver,level::HYPRE_Int)
    ccall((:HYPRE_GMRESSetPrintLevel,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,level)
end

function HYPRE_GMRESGetNumIterations(solver::HYPRE_Solver,num_iterations)
    ccall((:HYPRE_GMRESGetNumIterations,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,num_iterations)
end

function HYPRE_GMRESGetFinalRelativeResidualNorm(solver::HYPRE_Solver,norm)
    ccall((:HYPRE_GMRESGetFinalRelativeResidualNorm,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,norm)
end

function HYPRE_GMRESGetResidual(solver::HYPRE_Solver,residual)
    ccall((:HYPRE_GMRESGetResidual,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{Ptr{Void}}),solver,residual)
end

function HYPRE_GMRESGetTol(solver::HYPRE_Solver,tol)
    ccall((:HYPRE_GMRESGetTol,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,tol)
end

function HYPRE_GMRESGetAbsoluteTol(solver::HYPRE_Solver,tol)
    ccall((:HYPRE_GMRESGetAbsoluteTol,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,tol)
end

function HYPRE_GMRESGetConvergenceFactorTol(solver::HYPRE_Solver,cf_tol)
    ccall((:HYPRE_GMRESGetConvergenceFactorTol,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,cf_tol)
end

function HYPRE_GMRESGetStopCrit(solver::HYPRE_Solver,stop_crit)
    ccall((:HYPRE_GMRESGetStopCrit,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,stop_crit)
end

function HYPRE_GMRESGetMinIter(solver::HYPRE_Solver,min_iter)
    ccall((:HYPRE_GMRESGetMinIter,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,min_iter)
end

function HYPRE_GMRESGetMaxIter(solver::HYPRE_Solver,max_iter)
    ccall((:HYPRE_GMRESGetMaxIter,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,max_iter)
end

function HYPRE_GMRESGetKDim(solver::HYPRE_Solver,k_dim)
    ccall((:HYPRE_GMRESGetKDim,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,k_dim)
end

function HYPRE_GMRESGetRelChange(solver::HYPRE_Solver,rel_change)
    ccall((:HYPRE_GMRESGetRelChange,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,rel_change)
end

function HYPRE_GMRESGetPrecond(solver::HYPRE_Solver,precond_data_ptr)
    ccall((:HYPRE_GMRESGetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Solver}),solver,precond_data_ptr)
end

function HYPRE_GMRESGetLogging(solver::HYPRE_Solver,level)
    ccall((:HYPRE_GMRESGetLogging,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,level)
end

function HYPRE_GMRESGetPrintLevel(solver::HYPRE_Solver,level)
    ccall((:HYPRE_GMRESGetPrintLevel,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,level)
end

function HYPRE_GMRESGetConverged(solver::HYPRE_Solver,converged)
    ccall((:HYPRE_GMRESGetConverged,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,converged)
end

function HYPRE_FlexGMRESSetup(solver::HYPRE_Solver,A::HYPRE_Matrix,b::HYPRE_Vector,x::HYPRE_Vector)
    ccall((:HYPRE_FlexGMRESSetup,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector,HYPRE_Vector),solver,A,b,x)
end

function HYPRE_FlexGMRESSolve(solver::HYPRE_Solver,A::HYPRE_Matrix,b::HYPRE_Vector,x::HYPRE_Vector)
    ccall((:HYPRE_FlexGMRESSolve,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector,HYPRE_Vector),solver,A,b,x)
end

function HYPRE_FlexGMRESSetTol(solver::HYPRE_Solver,tol::HYPRE_Real)
    ccall((:HYPRE_FlexGMRESSetTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,tol)
end

function HYPRE_FlexGMRESSetAbsoluteTol(solver::HYPRE_Solver,a_tol::HYPRE_Real)
    ccall((:HYPRE_FlexGMRESSetAbsoluteTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,a_tol)
end

function HYPRE_FlexGMRESSetConvergenceFactorTol(solver::HYPRE_Solver,cf_tol::HYPRE_Real)
    ccall((:HYPRE_FlexGMRESSetConvergenceFactorTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,cf_tol)
end

function HYPRE_FlexGMRESSetMinIter(solver::HYPRE_Solver,min_iter::HYPRE_Int)
    ccall((:HYPRE_FlexGMRESSetMinIter,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,min_iter)
end

function HYPRE_FlexGMRESSetMaxIter(solver::HYPRE_Solver,max_iter::HYPRE_Int)
    ccall((:HYPRE_FlexGMRESSetMaxIter,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,max_iter)
end

function HYPRE_FlexGMRESSetKDim(solver::HYPRE_Solver,k_dim::HYPRE_Int)
    ccall((:HYPRE_FlexGMRESSetKDim,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,k_dim)
end

function HYPRE_FlexGMRESSetPrecond(solver::HYPRE_Solver,precond::HYPRE_PtrToSolverFcn,precond_setup::HYPRE_PtrToSolverFcn,precond_solver::HYPRE_Solver)
    ccall((:HYPRE_FlexGMRESSetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_PtrToSolverFcn,HYPRE_PtrToSolverFcn,HYPRE_Solver),solver,precond,precond_setup,precond_solver)
end

function HYPRE_FlexGMRESSetLogging(solver::HYPRE_Solver,logging::HYPRE_Int)
    ccall((:HYPRE_FlexGMRESSetLogging,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,logging)
end

function HYPRE_FlexGMRESSetPrintLevel(solver::HYPRE_Solver,level::HYPRE_Int)
    ccall((:HYPRE_FlexGMRESSetPrintLevel,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,level)
end

function HYPRE_FlexGMRESGetNumIterations(solver::HYPRE_Solver,num_iterations)
    ccall((:HYPRE_FlexGMRESGetNumIterations,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,num_iterations)
end

function HYPRE_FlexGMRESGetFinalRelativeResidualNorm(solver::HYPRE_Solver,norm)
    ccall((:HYPRE_FlexGMRESGetFinalRelativeResidualNorm,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,norm)
end

function HYPRE_FlexGMRESGetResidual(solver::HYPRE_Solver,residual)
    ccall((:HYPRE_FlexGMRESGetResidual,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{Ptr{Void}}),solver,residual)
end

function HYPRE_FlexGMRESGetTol(solver::HYPRE_Solver,tol)
    ccall((:HYPRE_FlexGMRESGetTol,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,tol)
end

function HYPRE_FlexGMRESGetConvergenceFactorTol(solver::HYPRE_Solver,cf_tol)
    ccall((:HYPRE_FlexGMRESGetConvergenceFactorTol,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,cf_tol)
end

function HYPRE_FlexGMRESGetStopCrit(solver::HYPRE_Solver,stop_crit)
    ccall((:HYPRE_FlexGMRESGetStopCrit,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,stop_crit)
end

function HYPRE_FlexGMRESGetMinIter(solver::HYPRE_Solver,min_iter)
    ccall((:HYPRE_FlexGMRESGetMinIter,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,min_iter)
end

function HYPRE_FlexGMRESGetMaxIter(solver::HYPRE_Solver,max_iter)
    ccall((:HYPRE_FlexGMRESGetMaxIter,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,max_iter)
end

function HYPRE_FlexGMRESGetKDim(solver::HYPRE_Solver,k_dim)
    ccall((:HYPRE_FlexGMRESGetKDim,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,k_dim)
end

function HYPRE_FlexGMRESGetPrecond(solver::HYPRE_Solver,precond_data_ptr)
    ccall((:HYPRE_FlexGMRESGetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Solver}),solver,precond_data_ptr)
end

function HYPRE_FlexGMRESGetLogging(solver::HYPRE_Solver,level)
    ccall((:HYPRE_FlexGMRESGetLogging,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,level)
end

function HYPRE_FlexGMRESGetPrintLevel(solver::HYPRE_Solver,level)
    ccall((:HYPRE_FlexGMRESGetPrintLevel,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,level)
end

function HYPRE_FlexGMRESGetConverged(solver::HYPRE_Solver,converged)
    ccall((:HYPRE_FlexGMRESGetConverged,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,converged)
end

function HYPRE_FlexGMRESSetModifyPC(solver::HYPRE_Solver,modify_pc::HYPRE_PtrToModifyPCFcn)
    ccall((:HYPRE_FlexGMRESSetModifyPC,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_PtrToModifyPCFcn),solver,modify_pc)
end

function HYPRE_LGMRESSetup(solver::HYPRE_Solver,A::HYPRE_Matrix,b::HYPRE_Vector,x::HYPRE_Vector)
    ccall((:HYPRE_LGMRESSetup,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector,HYPRE_Vector),solver,A,b,x)
end

function HYPRE_LGMRESSolve(solver::HYPRE_Solver,A::HYPRE_Matrix,b::HYPRE_Vector,x::HYPRE_Vector)
    ccall((:HYPRE_LGMRESSolve,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector,HYPRE_Vector),solver,A,b,x)
end

function HYPRE_LGMRESSetTol(solver::HYPRE_Solver,tol::HYPRE_Real)
    ccall((:HYPRE_LGMRESSetTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,tol)
end

function HYPRE_LGMRESSetAbsoluteTol(solver::HYPRE_Solver,a_tol::HYPRE_Real)
    ccall((:HYPRE_LGMRESSetAbsoluteTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,a_tol)
end

function HYPRE_LGMRESSetConvergenceFactorTol(solver::HYPRE_Solver,cf_tol::HYPRE_Real)
    ccall((:HYPRE_LGMRESSetConvergenceFactorTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,cf_tol)
end

function HYPRE_LGMRESSetMinIter(solver::HYPRE_Solver,min_iter::HYPRE_Int)
    ccall((:HYPRE_LGMRESSetMinIter,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,min_iter)
end

function HYPRE_LGMRESSetMaxIter(solver::HYPRE_Solver,max_iter::HYPRE_Int)
    ccall((:HYPRE_LGMRESSetMaxIter,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,max_iter)
end

function HYPRE_LGMRESSetKDim(solver::HYPRE_Solver,k_dim::HYPRE_Int)
    ccall((:HYPRE_LGMRESSetKDim,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,k_dim)
end

function HYPRE_LGMRESSetAugDim(solver::HYPRE_Solver,aug_dim::HYPRE_Int)
    ccall((:HYPRE_LGMRESSetAugDim,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,aug_dim)
end

function HYPRE_LGMRESSetPrecond(solver::HYPRE_Solver,precond::HYPRE_PtrToSolverFcn,precond_setup::HYPRE_PtrToSolverFcn,precond_solver::HYPRE_Solver)
    ccall((:HYPRE_LGMRESSetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_PtrToSolverFcn,HYPRE_PtrToSolverFcn,HYPRE_Solver),solver,precond,precond_setup,precond_solver)
end

function HYPRE_LGMRESSetLogging(solver::HYPRE_Solver,logging::HYPRE_Int)
    ccall((:HYPRE_LGMRESSetLogging,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,logging)
end

function HYPRE_LGMRESSetPrintLevel(solver::HYPRE_Solver,level::HYPRE_Int)
    ccall((:HYPRE_LGMRESSetPrintLevel,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,level)
end

function HYPRE_LGMRESGetNumIterations(solver::HYPRE_Solver,num_iterations)
    ccall((:HYPRE_LGMRESGetNumIterations,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,num_iterations)
end

function HYPRE_LGMRESGetFinalRelativeResidualNorm(solver::HYPRE_Solver,norm)
    ccall((:HYPRE_LGMRESGetFinalRelativeResidualNorm,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,norm)
end

function HYPRE_LGMRESGetResidual(solver::HYPRE_Solver,residual)
    ccall((:HYPRE_LGMRESGetResidual,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{Ptr{Void}}),solver,residual)
end

function HYPRE_LGMRESGetTol(solver::HYPRE_Solver,tol)
    ccall((:HYPRE_LGMRESGetTol,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,tol)
end

function HYPRE_LGMRESGetConvergenceFactorTol(solver::HYPRE_Solver,cf_tol)
    ccall((:HYPRE_LGMRESGetConvergenceFactorTol,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,cf_tol)
end

function HYPRE_LGMRESGetStopCrit(solver::HYPRE_Solver,stop_crit)
    ccall((:HYPRE_LGMRESGetStopCrit,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,stop_crit)
end

function HYPRE_LGMRESGetMinIter(solver::HYPRE_Solver,min_iter)
    ccall((:HYPRE_LGMRESGetMinIter,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,min_iter)
end

function HYPRE_LGMRESGetMaxIter(solver::HYPRE_Solver,max_iter)
    ccall((:HYPRE_LGMRESGetMaxIter,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,max_iter)
end

function HYPRE_LGMRESGetKDim(solver::HYPRE_Solver,k_dim)
    ccall((:HYPRE_LGMRESGetKDim,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,k_dim)
end

function HYPRE_LGMRESGetAugDim(solver::HYPRE_Solver,k_dim)
    ccall((:HYPRE_LGMRESGetAugDim,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,k_dim)
end

function HYPRE_LGMRESGetPrecond(solver::HYPRE_Solver,precond_data_ptr)
    ccall((:HYPRE_LGMRESGetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Solver}),solver,precond_data_ptr)
end

function HYPRE_LGMRESGetLogging(solver::HYPRE_Solver,level)
    ccall((:HYPRE_LGMRESGetLogging,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,level)
end

function HYPRE_LGMRESGetPrintLevel(solver::HYPRE_Solver,level)
    ccall((:HYPRE_LGMRESGetPrintLevel,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,level)
end

function HYPRE_LGMRESGetConverged(solver::HYPRE_Solver,converged)
    ccall((:HYPRE_LGMRESGetConverged,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,converged)
end

function HYPRE_BiCGSTABDestroy(solver::HYPRE_Solver)
    ccall((:HYPRE_BiCGSTABDestroy,libHypre),HYPRE_Int,(HYPRE_Solver,),solver)
end

function HYPRE_BiCGSTABSetup(solver::HYPRE_Solver,A::HYPRE_Matrix,b::HYPRE_Vector,x::HYPRE_Vector)
    ccall((:HYPRE_BiCGSTABSetup,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector,HYPRE_Vector),solver,A,b,x)
end

function HYPRE_BiCGSTABSolve(solver::HYPRE_Solver,A::HYPRE_Matrix,b::HYPRE_Vector,x::HYPRE_Vector)
    ccall((:HYPRE_BiCGSTABSolve,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector,HYPRE_Vector),solver,A,b,x)
end

function HYPRE_BiCGSTABSetTol(solver::HYPRE_Solver,tol::HYPRE_Real)
    ccall((:HYPRE_BiCGSTABSetTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,tol)
end

function HYPRE_BiCGSTABSetAbsoluteTol(solver::HYPRE_Solver,a_tol::HYPRE_Real)
    ccall((:HYPRE_BiCGSTABSetAbsoluteTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,a_tol)
end

function HYPRE_BiCGSTABSetConvergenceFactorTol(solver::HYPRE_Solver,cf_tol::HYPRE_Real)
    ccall((:HYPRE_BiCGSTABSetConvergenceFactorTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,cf_tol)
end

function HYPRE_BiCGSTABSetStopCrit(solver::HYPRE_Solver,stop_crit::HYPRE_Int)
    ccall((:HYPRE_BiCGSTABSetStopCrit,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,stop_crit)
end

function HYPRE_BiCGSTABSetMinIter(solver::HYPRE_Solver,min_iter::HYPRE_Int)
    ccall((:HYPRE_BiCGSTABSetMinIter,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,min_iter)
end

function HYPRE_BiCGSTABSetMaxIter(solver::HYPRE_Solver,max_iter::HYPRE_Int)
    ccall((:HYPRE_BiCGSTABSetMaxIter,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,max_iter)
end

function HYPRE_BiCGSTABSetPrecond(solver::HYPRE_Solver,precond::HYPRE_PtrToSolverFcn,precond_setup::HYPRE_PtrToSolverFcn,precond_solver::HYPRE_Solver)
    ccall((:HYPRE_BiCGSTABSetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_PtrToSolverFcn,HYPRE_PtrToSolverFcn,HYPRE_Solver),solver,precond,precond_setup,precond_solver)
end

function HYPRE_BiCGSTABSetLogging(solver::HYPRE_Solver,logging::HYPRE_Int)
    ccall((:HYPRE_BiCGSTABSetLogging,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,logging)
end

function HYPRE_BiCGSTABSetPrintLevel(solver::HYPRE_Solver,level::HYPRE_Int)
    ccall((:HYPRE_BiCGSTABSetPrintLevel,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,level)
end

function HYPRE_BiCGSTABGetNumIterations(solver::HYPRE_Solver,num_iterations)
    ccall((:HYPRE_BiCGSTABGetNumIterations,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,num_iterations)
end

function HYPRE_BiCGSTABGetFinalRelativeResidualNorm(solver::HYPRE_Solver,norm)
    ccall((:HYPRE_BiCGSTABGetFinalRelativeResidualNorm,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,norm)
end

function HYPRE_BiCGSTABGetResidual(solver::HYPRE_Solver,residual)
    ccall((:HYPRE_BiCGSTABGetResidual,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{Ptr{Void}}),solver,residual)
end

function HYPRE_BiCGSTABGetPrecond(solver::HYPRE_Solver,precond_data_ptr)
    ccall((:HYPRE_BiCGSTABGetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Solver}),solver,precond_data_ptr)
end

function HYPRE_CGNRDestroy(solver::HYPRE_Solver)
    ccall((:HYPRE_CGNRDestroy,libHypre),HYPRE_Int,(HYPRE_Solver,),solver)
end

function HYPRE_CGNRSetup(solver::HYPRE_Solver,A::HYPRE_Matrix,b::HYPRE_Vector,x::HYPRE_Vector)
    ccall((:HYPRE_CGNRSetup,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector,HYPRE_Vector),solver,A,b,x)
end

function HYPRE_CGNRSolve(solver::HYPRE_Solver,A::HYPRE_Matrix,b::HYPRE_Vector,x::HYPRE_Vector)
    ccall((:HYPRE_CGNRSolve,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector,HYPRE_Vector),solver,A,b,x)
end

function HYPRE_CGNRSetTol(solver::HYPRE_Solver,tol::HYPRE_Real)
    ccall((:HYPRE_CGNRSetTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,tol)
end

function HYPRE_CGNRSetStopCrit(solver::HYPRE_Solver,stop_crit::HYPRE_Int)
    ccall((:HYPRE_CGNRSetStopCrit,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,stop_crit)
end

function HYPRE_CGNRSetMinIter(solver::HYPRE_Solver,min_iter::HYPRE_Int)
    ccall((:HYPRE_CGNRSetMinIter,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,min_iter)
end

function HYPRE_CGNRSetMaxIter(solver::HYPRE_Solver,max_iter::HYPRE_Int)
    ccall((:HYPRE_CGNRSetMaxIter,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,max_iter)
end

function HYPRE_CGNRSetPrecond(solver::HYPRE_Solver,precond::HYPRE_PtrToSolverFcn,precondT::HYPRE_PtrToSolverFcn,precond_setup::HYPRE_PtrToSolverFcn,precond_solver::HYPRE_Solver)
    ccall((:HYPRE_CGNRSetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_PtrToSolverFcn,HYPRE_PtrToSolverFcn,HYPRE_PtrToSolverFcn,HYPRE_Solver),solver,precond,precondT,precond_setup,precond_solver)
end

function HYPRE_CGNRSetLogging(solver::HYPRE_Solver,logging::HYPRE_Int)
    ccall((:HYPRE_CGNRSetLogging,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,logging)
end

function HYPRE_CGNRGetNumIterations(solver::HYPRE_Solver,num_iterations)
    ccall((:HYPRE_CGNRGetNumIterations,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Int}),solver,num_iterations)
end

function HYPRE_CGNRGetFinalRelativeResidualNorm(solver::HYPRE_Solver,norm)
    ccall((:HYPRE_CGNRGetFinalRelativeResidualNorm,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Real}),solver,norm)
end

function HYPRE_CGNRGetPrecond(solver::HYPRE_Solver,precond_data_ptr)
    ccall((:HYPRE_CGNRGetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Solver}),solver,precond_data_ptr)
end

function utilities_FortranMatrixCreate()
    ccall((:utilities_FortranMatrixCreate,libHypre),Ptr{utilities_FortranMatrix},())
end

function utilities_FortranMatrixAllocateData(h::hypre_longint,w::hypre_longint,mtx)
    ccall((:utilities_FortranMatrixAllocateData,libHypre),Void,(hypre_longint,hypre_longint,Ptr{utilities_FortranMatrix}),h,w,mtx)
end

function utilities_FortranMatrixWrap(arg1,gh::hypre_longint,h::hypre_longint,w::hypre_longint,mtx)
    ccall((:utilities_FortranMatrixWrap,libHypre),Void,(Ptr{HYPRE_Real},hypre_longint,hypre_longint,hypre_longint,Ptr{utilities_FortranMatrix}),arg1,gh,h,w,mtx)
end

function utilities_FortranMatrixDestroy(mtx)
    ccall((:utilities_FortranMatrixDestroy,libHypre),Void,(Ptr{utilities_FortranMatrix},),mtx)
end

function utilities_FortranMatrixGlobalHeight(mtx)
    ccall((:utilities_FortranMatrixGlobalHeight,libHypre),hypre_longint,(Ptr{utilities_FortranMatrix},),mtx)
end

function utilities_FortranMatrixHeight(mtx)
    ccall((:utilities_FortranMatrixHeight,libHypre),hypre_longint,(Ptr{utilities_FortranMatrix},),mtx)
end

function utilities_FortranMatrixWidth(mtx)
    ccall((:utilities_FortranMatrixWidth,libHypre),hypre_longint,(Ptr{utilities_FortranMatrix},),mtx)
end

function utilities_FortranMatrixValues(mtx)
    ccall((:utilities_FortranMatrixValues,libHypre),Ptr{HYPRE_Real},(Ptr{utilities_FortranMatrix},),mtx)
end

function utilities_FortranMatrixClear(mtx)
    ccall((:utilities_FortranMatrixClear,libHypre),Void,(Ptr{utilities_FortranMatrix},),mtx)
end

function utilities_FortranMatrixClearL(mtx)
    ccall((:utilities_FortranMatrixClearL,libHypre),Void,(Ptr{utilities_FortranMatrix},),mtx)
end

function utilities_FortranMatrixSetToIdentity(mtx)
    ccall((:utilities_FortranMatrixSetToIdentity,libHypre),Void,(Ptr{utilities_FortranMatrix},),mtx)
end

function utilities_FortranMatrixTransposeSquare(mtx)
    ccall((:utilities_FortranMatrixTransposeSquare,libHypre),Void,(Ptr{utilities_FortranMatrix},),mtx)
end

function utilities_FortranMatrixSymmetrize(mtx)
    ccall((:utilities_FortranMatrixSymmetrize,libHypre),Void,(Ptr{utilities_FortranMatrix},),mtx)
end

function utilities_FortranMatrixCopy(src,t::HYPRE_Int,dest)
    ccall((:utilities_FortranMatrixCopy,libHypre),Void,(Ptr{utilities_FortranMatrix},HYPRE_Int,Ptr{utilities_FortranMatrix}),src,t,dest)
end

function utilities_FortranMatrixIndexCopy(index,src,t::HYPRE_Int,dest)
    ccall((:utilities_FortranMatrixIndexCopy,libHypre),Void,(Ptr{HYPRE_Int},Ptr{utilities_FortranMatrix},HYPRE_Int,Ptr{utilities_FortranMatrix}),index,src,t,dest)
end

function utilities_FortranMatrixSetDiagonal(mtx,d)
    ccall((:utilities_FortranMatrixSetDiagonal,libHypre),Void,(Ptr{utilities_FortranMatrix},Ptr{utilities_FortranMatrix}),mtx,d)
end

function utilities_FortranMatrixGetDiagonal(mtx,d)
    ccall((:utilities_FortranMatrixGetDiagonal,libHypre),Void,(Ptr{utilities_FortranMatrix},Ptr{utilities_FortranMatrix}),mtx,d)
end

function utilities_FortranMatrixAdd(a::HYPRE_Real,mtxA,mtxB,mtxC)
    ccall((:utilities_FortranMatrixAdd,libHypre),Void,(HYPRE_Real,Ptr{utilities_FortranMatrix},Ptr{utilities_FortranMatrix},Ptr{utilities_FortranMatrix}),a,mtxA,mtxB,mtxC)
end

function utilities_FortranMatrixDMultiply(d,mtx)
    ccall((:utilities_FortranMatrixDMultiply,libHypre),Void,(Ptr{utilities_FortranMatrix},Ptr{utilities_FortranMatrix}),d,mtx)
end

function utilities_FortranMatrixMultiplyD(mtx,d)
    ccall((:utilities_FortranMatrixMultiplyD,libHypre),Void,(Ptr{utilities_FortranMatrix},Ptr{utilities_FortranMatrix}),mtx,d)
end

function utilities_FortranMatrixMultiply(mtxA,tA::HYPRE_Int,mtxB,tB::HYPRE_Int,mtxC)
    ccall((:utilities_FortranMatrixMultiply,libHypre),Void,(Ptr{utilities_FortranMatrix},HYPRE_Int,Ptr{utilities_FortranMatrix},HYPRE_Int,Ptr{utilities_FortranMatrix}),mtxA,tA,mtxB,tB,mtxC)
end

function utilities_FortranMatrixFNorm(mtx)
    ccall((:utilities_FortranMatrixFNorm,libHypre),HYPRE_Real,(Ptr{utilities_FortranMatrix},),mtx)
end

function utilities_FortranMatrixValue(mtx,i::hypre_longint,j::hypre_longint)
    ccall((:utilities_FortranMatrixValue,libHypre),HYPRE_Real,(Ptr{utilities_FortranMatrix},hypre_longint,hypre_longint),mtx,i,j)
end

function utilities_FortranMatrixValuePtr(mtx,i::hypre_longint,j::hypre_longint)
    ccall((:utilities_FortranMatrixValuePtr,libHypre),Ptr{HYPRE_Real},(Ptr{utilities_FortranMatrix},hypre_longint,hypre_longint),mtx,i,j)
end

function utilities_FortranMatrixMaxValue(mtx)
    ccall((:utilities_FortranMatrixMaxValue,libHypre),HYPRE_Real,(Ptr{utilities_FortranMatrix},),mtx)
end

function utilities_FortranMatrixSelectBlock(mtx,iFrom::hypre_longint,iTo::hypre_longint,jFrom::hypre_longint,jTo::hypre_longint,block)
    ccall((:utilities_FortranMatrixSelectBlock,libHypre),Void,(Ptr{utilities_FortranMatrix},hypre_longint,hypre_longint,hypre_longint,hypre_longint,Ptr{utilities_FortranMatrix}),mtx,iFrom,iTo,jFrom,jTo,block)
end

function utilities_FortranMatrixUpperInv(u)
    ccall((:utilities_FortranMatrixUpperInv,libHypre),Void,(Ptr{utilities_FortranMatrix},),u)
end

function utilities_FortranMatrixPrint(mtx,fileName)
    ccall((:utilities_FortranMatrixPrint,libHypre),HYPRE_Int,(Ptr{utilities_FortranMatrix},Cstring),mtx,fileName)
end

function mv_MultiVectorGetData(x::mv_MultiVectorPtr)
    ccall((:mv_MultiVectorGetData,libHypre),Ptr{Void},(mv_MultiVectorPtr,),x)
end

function mv_MultiVectorWrap(ii,data,ownsData::HYPRE_Int)
    ccall((:mv_MultiVectorWrap,libHypre),mv_MultiVectorPtr,(Ptr{mv_InterfaceInterpreter},Ptr{Void},HYPRE_Int),ii,data,ownsData)
end

function mv_MultiVectorCreateFromSampleVector(arg1,n::HYPRE_Int,sample)
    ccall((:mv_MultiVectorCreateFromSampleVector,libHypre),mv_MultiVectorPtr,(Ptr{Void},HYPRE_Int,Ptr{Void}),arg1,n,sample)
end

function mv_MultiVectorCreateCopy(x::mv_MultiVectorPtr,copyValues::HYPRE_Int)
    ccall((:mv_MultiVectorCreateCopy,libHypre),mv_MultiVectorPtr,(mv_MultiVectorPtr,HYPRE_Int),x,copyValues)
end

function mv_MultiVectorDestroy(arg1::mv_MultiVectorPtr)
    ccall((:mv_MultiVectorDestroy,libHypre),Void,(mv_MultiVectorPtr,),arg1)
end

function mv_MultiVectorWidth(v::mv_MultiVectorPtr)
    ccall((:mv_MultiVectorWidth,libHypre),HYPRE_Int,(mv_MultiVectorPtr,),v)
end

function mv_MultiVectorHeight(v::mv_MultiVectorPtr)
    ccall((:mv_MultiVectorHeight,libHypre),HYPRE_Int,(mv_MultiVectorPtr,),v)
end

function mv_MultiVectorSetMask(v::mv_MultiVectorPtr,mask)
    ccall((:mv_MultiVectorSetMask,libHypre),Void,(mv_MultiVectorPtr,Ptr{HYPRE_Int}),v,mask)
end

function mv_MultiVectorClear(arg1::mv_MultiVectorPtr)
    ccall((:mv_MultiVectorClear,libHypre),Void,(mv_MultiVectorPtr,),arg1)
end

function mv_MultiVectorSetRandom(v::mv_MultiVectorPtr,seed::HYPRE_Int)
    ccall((:mv_MultiVectorSetRandom,libHypre),Void,(mv_MultiVectorPtr,HYPRE_Int),v,seed)
end

function mv_MultiVectorCopy(src::mv_MultiVectorPtr,dest::mv_MultiVectorPtr)
    ccall((:mv_MultiVectorCopy,libHypre),Void,(mv_MultiVectorPtr,mv_MultiVectorPtr),src,dest)
end

function mv_MultiVectorAxpy(a::HYPRE_Complex,x::mv_MultiVectorPtr,y::mv_MultiVectorPtr)
    ccall((:mv_MultiVectorAxpy,libHypre),Void,(HYPRE_Complex,mv_MultiVectorPtr,mv_MultiVectorPtr),a,x,y)
end

function mv_MultiVectorByMultiVector(x::mv_MultiVectorPtr,y::mv_MultiVectorPtr,gh::HYPRE_Int,h::HYPRE_Int,w::HYPRE_Int,v)
    ccall((:mv_MultiVectorByMultiVector,libHypre),Void,(mv_MultiVectorPtr,mv_MultiVectorPtr,HYPRE_Int,HYPRE_Int,HYPRE_Int,Ptr{HYPRE_Real}),x,y,gh,h,w,v)
end

function mv_MultiVectorByMultiVectorDiag(arg1::mv_MultiVectorPtr,arg2::mv_MultiVectorPtr,mask,n::HYPRE_Int,diag)
    ccall((:mv_MultiVectorByMultiVectorDiag,libHypre),Void,(mv_MultiVectorPtr,mv_MultiVectorPtr,Ptr{HYPRE_Int},HYPRE_Int,Ptr{HYPRE_Real}),arg1,arg2,mask,n,diag)
end

function mv_MultiVectorByMatrix(x::mv_MultiVectorPtr,gh::HYPRE_Int,h::HYPRE_Int,w::HYPRE_Int,v,y::mv_MultiVectorPtr)
    ccall((:mv_MultiVectorByMatrix,libHypre),Void,(mv_MultiVectorPtr,HYPRE_Int,HYPRE_Int,HYPRE_Int,Ptr{HYPRE_Complex},mv_MultiVectorPtr),x,gh,h,w,v,y)
end

function mv_MultiVectorXapy(x::mv_MultiVectorPtr,gh::HYPRE_Int,h::HYPRE_Int,w::HYPRE_Int,v,y::mv_MultiVectorPtr)
    ccall((:mv_MultiVectorXapy,libHypre),Void,(mv_MultiVectorPtr,HYPRE_Int,HYPRE_Int,HYPRE_Int,Ptr{HYPRE_Complex},mv_MultiVectorPtr),x,gh,h,w,v,y)
end

function mv_MultiVectorByDiagonal(x::mv_MultiVectorPtr,mask,n::HYPRE_Int,diag,y::mv_MultiVectorPtr)
    ccall((:mv_MultiVectorByDiagonal,libHypre),Void,(mv_MultiVectorPtr,Ptr{HYPRE_Int},HYPRE_Int,Ptr{HYPRE_Complex},mv_MultiVectorPtr),x,mask,n,diag,y)
end

function mv_MultiVectorEval(f,par,x::mv_MultiVectorPtr,y::mv_MultiVectorPtr)
    ccall((:mv_MultiVectorEval,libHypre),Void,(Ptr{Void},Ptr{Void},mv_MultiVectorPtr,mv_MultiVectorPtr),f,par,x,y)
end

function HYPRE_LOBPCGCreate(interpreter,mvfunctions,solver)
    ccall((:HYPRE_LOBPCGCreate,libHypre),HYPRE_Int,(Ptr{mv_InterfaceInterpreter},Ptr{HYPRE_MatvecFunctions},Ptr{HYPRE_Solver}),interpreter,mvfunctions,solver)
end

function HYPRE_LOBPCGDestroy(solver::HYPRE_Solver)
    ccall((:HYPRE_LOBPCGDestroy,libHypre),HYPRE_Int,(HYPRE_Solver,),solver)
end

function HYPRE_LOBPCGSetPrecond(solver::HYPRE_Solver,precond::HYPRE_PtrToSolverFcn,precond_setup::HYPRE_PtrToSolverFcn,precond_solver::HYPRE_Solver)
    ccall((:HYPRE_LOBPCGSetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_PtrToSolverFcn,HYPRE_PtrToSolverFcn,HYPRE_Solver),solver,precond,precond_setup,precond_solver)
end

function HYPRE_LOBPCGGetPrecond(solver::HYPRE_Solver,precond_data_ptr)
    ccall((:HYPRE_LOBPCGGetPrecond,libHypre),HYPRE_Int,(HYPRE_Solver,Ptr{HYPRE_Solver}),solver,precond_data_ptr)
end

function HYPRE_LOBPCGSetup(solver::HYPRE_Solver,A::HYPRE_Matrix,b::HYPRE_Vector,x::HYPRE_Vector)
    ccall((:HYPRE_LOBPCGSetup,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector,HYPRE_Vector),solver,A,b,x)
end

function HYPRE_LOBPCGSetupB(solver::HYPRE_Solver,B::HYPRE_Matrix,x::HYPRE_Vector)
    ccall((:HYPRE_LOBPCGSetupB,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector),solver,B,x)
end

function HYPRE_LOBPCGSetupT(solver::HYPRE_Solver,T::HYPRE_Matrix,x::HYPRE_Vector)
    ccall((:HYPRE_LOBPCGSetupT,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Matrix,HYPRE_Vector),solver,T,x)
end

function HYPRE_LOBPCGSolve(solver::HYPRE_Solver,y::mv_MultiVectorPtr,x::mv_MultiVectorPtr,lambda)
    ccall((:HYPRE_LOBPCGSolve,libHypre),HYPRE_Int,(HYPRE_Solver,mv_MultiVectorPtr,mv_MultiVectorPtr,Ptr{HYPRE_Real}),solver,y,x,lambda)
end

function HYPRE_LOBPCGSetTol(solver::HYPRE_Solver,tol::HYPRE_Real)
    ccall((:HYPRE_LOBPCGSetTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,tol)
end

function HYPRE_LOBPCGSetRTol(solver::HYPRE_Solver,tol::HYPRE_Real)
    ccall((:HYPRE_LOBPCGSetRTol,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Real),solver,tol)
end

function HYPRE_LOBPCGSetMaxIter(solver::HYPRE_Solver,max_iter::HYPRE_Int)
    ccall((:HYPRE_LOBPCGSetMaxIter,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,max_iter)
end

function HYPRE_LOBPCGSetPrecondUsageMode(solver::HYPRE_Solver,mode::HYPRE_Int)
    ccall((:HYPRE_LOBPCGSetPrecondUsageMode,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,mode)
end

function HYPRE_LOBPCGSetPrintLevel(solver::HYPRE_Solver,level::HYPRE_Int)
    ccall((:HYPRE_LOBPCGSetPrintLevel,libHypre),HYPRE_Int,(HYPRE_Solver,HYPRE_Int),solver,level)
end

function HYPRE_LOBPCGResidualNorms(solver::HYPRE_Solver)
    ccall((:HYPRE_LOBPCGResidualNorms,libHypre),Ptr{utilities_FortranMatrix},(HYPRE_Solver,),solver)
end

function HYPRE_LOBPCGResidualNormsHistory(solver::HYPRE_Solver)
    ccall((:HYPRE_LOBPCGResidualNormsHistory,libHypre),Ptr{utilities_FortranMatrix},(HYPRE_Solver,),solver)
end

function HYPRE_LOBPCGEigenvaluesHistory(solver::HYPRE_Solver)
    ccall((:HYPRE_LOBPCGEigenvaluesHistory,libHypre),Ptr{utilities_FortranMatrix},(HYPRE_Solver,),solver)
end

function HYPRE_LOBPCGIterations(solver::HYPRE_Solver)
    ccall((:HYPRE_LOBPCGIterations,libHypre),HYPRE_Int,(HYPRE_Solver,),solver)
end

function hypre_LOBPCGMultiOperatorB(data,x,y)
    ccall((:hypre_LOBPCGMultiOperatorB,libHypre),Void,(Ptr{Void},Ptr{Void},Ptr{Void}),data,x,y)
end

function lobpcg_MultiVectorByMultiVector(x::mv_MultiVectorPtr,y::mv_MultiVectorPtr,xy)
    ccall((:lobpcg_MultiVectorByMultiVector,libHypre),Void,(mv_MultiVectorPtr,mv_MultiVectorPtr,Ptr{utilities_FortranMatrix}),x,y,xy)
end


text:
@testset "k_cliques tests" begin

    @testset "triangle tests" begin
        @test length(kcliques(MatrixNetwork(sparse([0 1; 1 0])),3)[2]) == 0
        @test length(kcliques(MatrixNetwork(sparse(fill(0,10,10))),3)[2]) == 0
        @test length(kcliques(MatrixNetwork(sparse(fill(0,10,10))),3)[2]) == 0
        @test length(kcliques(MatrixNetwork(sparse(fill(0,10,10))),3)[2]) == 0
        @test length(kcliques(MatrixNetwork(sparse([0 1 1; 1 0 1; 1 1 0])), 3)[2]) == 1
        @test length(kcliques(load_data("facebook_data"), 3)[2]) == 1612010
        @test length(kcliques(MatrixNetwork(load_matrix_network("clique-10")), 3)[2]) == 120
    end

end


text:
function _local_baseline(xs, ys, xₗ, xᵣ, b::UncertainBound)
    ŷₗ = mean(lininterp(q, xs, ys) for q in b._left_quantiles)
    ŷᵣ = mean(lininterp(q, xs, ys) for q in b._right_quantiles)
    x̂ₗ = b._left_quantiles[IDX_BOUND_MEDIAN]
    x̂ᵣ = b._right_quantiles[IDX_BOUND_MEDIAN]
    yₗ = lininterp(xₗ, x̂ₗ, x̂ᵣ, ŷₗ, ŷᵣ)
    yᵣ = lininterp(xᵣ, x̂ₗ, x̂ᵣ, ŷₗ, ŷᵣ)
    return xₗ, xᵣ, yₗ, yᵣ
end

_endpoint_to_endpoint_baseline(xs, ys, xₗ, xᵣ) = (xₗ, xᵣ, lininterp(xₗ, xs, ys), lininterp(xᵣ, xs, ys))

"""
    mc_integrate(uc::UncertainCurve{T, N}, bnds::Vector{UncertainBound{T, M}}; intfun=trapz)
    mc_integrate(uc::UncertainCurve{T, N}, bnds::UncertainBound{T, M}; intfun=trapz)

Integrate uncertain curve using uncertain bound(s).

Applies Monte-Carlo integration algorithm to samples stored in `uc` and returns one or an array
of uncertain areas of type `Particles{T, N}`.

# Keyword arguments

`intfun`:
The core integration function that is used to numerically integrate each draw.
Defaults to `NoisySignalIntegration.trapz`.
The function that is used to substitute [`trapz`](@ref) must share its call signature.

`subtract_baseline` (deprecated in favor of `local_baseline`):
If true, for each draw a local linear baseline defined by the integration window start and end point will be subtracted.

`local_baseline`:
If true, for each draw a local linear baseline defined by the integration window start and end point will be subtracted.
The y-values of the start and end point are derived from a weighted average over the start and end point distributions, see
[the documentation](https://nluetts.github.io/NoisySignalIntegration.jl/dev/baseline/#Build-in) for further information.
"""
function mc_integrate(uc::UncertainCurve{T, N}, bnds::Vector{UncertainBound{T, M}}; intfun=trapz, subtract_baseline=false, local_baseline=false) where {T, M, N}

    M != N && error("Samples sizes incompatible")
    subtract_baseline && @warn("subtract_baseline keyword argument is deprecated, use local_baseline instead.")
    (subtract_baseline && local_baseline) && error("local_baseline and subtract_baseline cannot both be true.") |> throw

    areas = Array{T}(undef, N, length(bnds))
    for i ∈ 1:N
        i % 1000 == 0 && print("Integrating draw $i/$N \r")
        cᵢ = get_draw(i, uc)
        for (j, b) in enumerate(bnds)
            xₗ, xᵣ = get_draw(i, b)
            xs, ys = uc.x, cᵢ.y
            areas[i, j] = intfun(xs, ys, xₗ, xᵣ)
            if local_baseline
                areas[i, j] -= singletrapz(_local_baseline(xs, ys, xₗ, xᵣ, b)...)
            end
            if subtract_baseline
                areas[i, j] -= singletrapz(_endpoint_to_endpoint_baseline(xs, ys, xₗ, xᵣ)...)
            end
        end
    end
    return [Particles(areas[:,i]) for i in 1:size(areas)[2]]
end

function mc_integrate(uc::S, bnd::T; intfun=trapz, subtract_baseline=false, local_baseline=false) where {S <: UncertainCurve, T <: UncertainBound}
    return mc_integrate(uc, [bnd]; intfun=intfun, subtract_baseline=subtract_baseline, local_baseline=local_baseline)[1]
end


text:
# peering into InSight data
using CSV
using DataFrames
function read_temp_wind_data()
    # data found at https://atmos.nmsu.edu/data_and_services/atmospheres_data/INSIGHT/insight.html#Selecting_Data
    # placed in parent dir/data
    path = "../data/twins_bundle/data_calibrated/sol_0000_0122/"
    x  = readdir(path)
    fn = "twins_calibevent_0118_01.csv"
    df = CSV.read(string(path, fn))
    return df
end


text:
s = "hello"
println(s * " there!")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using DiffEqBase

####################################################################################################
# Poincare shooting based on Sánchez, J., M. Net, B. Garcı́a-Archilla, and C. Simó. “Newton–Krylov Continuation of Periodic Orbits for Navier–Stokes Flows.” Journal of Computational Physics 201, no. 1 (November 20, 2004): 13–33. https://doi.org/10.1016/j.jcp.2004.04.018.

function sectionHyp!(out, x, normals, centers)
	for ii = 1:length(normals)
		out[ii] = dot(normals[ii], x) - dot(normals[ii], centers[ii])
	end
	out
end

# this composite type encodes a set of hyperplanes which are used as Poincaré sections
struct HyperplaneSections{Tn, Tc, Ti, Tnb, Tcb}
	M::Int64		# number of hyperplanes
	normals::Tn 	# normals to define hyperplanes
	centers::Tc 	# representative point on each hyperplane
	indices::Ti 	# indices to be removed in the operator Ek

	normals_bar::Tnb
	centers_bar::Tcb

	function HyperplaneSections(normals, centers)
		M = length(normals)
		indices = zeros(Int64, M)
		for ii in 1:M
			indices[ii] = argmax(abs.(normals[ii]))
		end
		nbar = [R(normals[ii], indices[ii]) for ii=1:M]
		cbar = [R(centers[ii], indices[ii]) for ii=1:M]

		return new{typeof(normals), typeof(centers), typeof(indices), typeof(nbar), typeof(cbar)}(M, normals, centers, indices, nbar, cbar)
	end

	HyperplaneSections(M = 0) = new{Nothing, Nothing, Nothing, Nothing, Nothing}(M, nothing, nothing, nothing, nothing, nothing)
end

(hyp::HyperplaneSections)(out, u) = sectionHyp!(out, u, hyp.normals, hyp.centers)

"""
	update!(hyp::HyperplaneSections, normals, centers)

Update the hyperplanes saved in `hyp`.
"""
function update!(hyp::HyperplaneSections, normals, centers)
	M = hyp.M
	@assert length(normals) == M "Wrong number of normals"
	@assert length(centers) == M "Wrong number of centers"
	hyp.normals .= normals
	hyp.centers .= centers
	for ii in 1:M
		hyp.indices[ii] = argmax(abs.(normals[ii]))
		R!(hyp.normals_bar[ii], normals[ii], hyp.indices[ii])
		R!(hyp.centers_bar[ii], centers[ii], hyp.indices[ii])
	end
	return hyp
end

# Operateur Rk from the paper above
function R!(out, x::AbstractVector, k::Int)
	@views out[1:k-1] .= x[1:k-1]
	@views out[k:end] .= x[k+1:end]
	return out
end

R!(hyp::HyperplaneSections, out, x::AbstractVector, ii::Int) = R!(out, x, hyp.indices[ii])

function R(x::AbstractVector, k::Int)
	out = similar(x, length(x) - 1)
	R!(out, x, k)
end

function R(hyp::HyperplaneSections, x::AbstractVector, ii::Int)
	out = similar(x, length(x) - 1)
	R!(hyp, out, x, ii)
end

# differential of R
dR!(hyp::HyperplaneSections, out, dx::AbstractVector, ii::Int) = R!(hyp, out, dx, ii)

# Operateur Ek from the paper above
function E!(hyp::HyperplaneSections, out, xbar::AbstractVector, ii::Int)
	@assert length(xbar) == length(hyp.normals[1]) - 1 "Wrong size for the projector / expansion operators, length(xbar) = $(length(xbar)) and length(normal) = $(length(hyp.normals[1]))"
	k = hyp.indices[ii]
	nbar = hyp.normals_bar[ii]
	xcbar = hyp.centers_bar[ii]
	coord_k = hyp.centers[ii][k] - (dot(nbar, xbar) - dot(nbar, xcbar)) / hyp.normals[ii][k]

	@views out[1:k-1] .= xbar[1:k-1]
	@views out[k+1:end] .= xbar[k:end]
	out[k] = coord_k
	return out
end

function E(hyp::HyperplaneSections, xbar::AbstractVector, ii::Int)
	out = similar(xbar, length(xbar) + 1)
	E!(hyp, out, xbar, ii)
end

# differential of E!
function dE!(hyp::HyperplaneSections, out, dxbar::AbstractVector, ii::Int)
	k = hyp.indices[ii]
	nbar = hyp.normals_bar[ii]
	xcbar = hyp.centers_bar[ii]
	coord_k = - dot(nbar, dxbar) / hyp.normals[ii][k]

	@views out[1:k-1]   .= dxbar[1:k-1]
	@views out[k+1:end] .= dxbar[k:end]
	out[k] = coord_k
	return out
end

function dE(hyp::HyperplaneSections, dxbar::AbstractVector, ii::Int)
	out = similar(dxbar, length(dxbar) + 1)
	dE!(hyp, out, dxbar, ii)
end
####################################################################################################
"""

`pb = PoincareShootingProblem(flow::Flow, M, sections; δ = 1e-8, interp_points = 50, isparallel = false)`

This composite type implements the Poincaré Shooting method to locate periodic orbits by relying on Poincaré return maps. The arguments are as follows
- `flow::Flow`: implements the flow of the Cauchy problem though the structure [`Flow`](@ref).
- `M`: the number of Poincaré sections. If `M==1`, then the simple shooting is implemented and the multiple one otherwise.
- `sections`: function or callable strict which implements a Poincaré section condition. The evaluation `sections(x)` must return a scalar number when `M==1`. Otherwise, one must implement a function `section(out, x)` which populates `out` with the `M` sections.
- `δ = 1e-8` used to compute the jacobian of the fonctional by finite differences. If set to `0`, an analytical expression of the jacobian is used instead.
- `interp_points = 50` number of interpolation point used to define the callback (to compute the hitting of the hyperplan section)
- `isparallel = false` whether the shooting are computed in parallel (threading). Only available through the use of Flows defined by `EnsembleProblem`.

## Simplified constructors
- A simpler way is to create a functional is

`pb = PoincareShootingProblem(F, p, prob::ODEProblem, alg, section; kwargs...)`

for simple shooting or

`pb = PoincareShootingProblem(F, p, prob::Union{ODEProblem, EnsembleProblem}, alg, M::Int, section; kwargs...)`

for multiple shooting . Here `F(x,p)` is the vector field, `p` is a parameter (to be passed to the vector and the flow), `prob` is an `Union{ODEProblem, EnsembleProblem}` which is used to create a flow using the ODE solver `alg` (for example `Tsit5()`). Finally, the arguments `kwargs` are passed to the ODE solver defining the flow. We refere to `DifferentialEquations.jl` for more information.

- Another convenient call is

`pb = PoincareShootingProblem(F, p, prob::Union{ODEProblem, EnsembleProblem}, alg, normals::AbstractVector, centers::AbstractVector; δ = 1e-8, kwargs...)`

where `normals` (resp. `centers`) is a list of normals (resp. centers) which defines a list of hyperplanes ``\\Sigma_i``. These hyperplanes are used to define partial Poincaré return maps.

## Computing the functionals
A functional, hereby called `G` encodes this shooting problem. You can then call `pb(orbitguess, par)` to apply the functional to a guess. Note that `orbitguess::AbstractVector` must be of size M * N where N is the number of unknowns in the state space and `M` is the number of Poincaré maps. Another accepted `guess` is such that `guess[i]` is the state of the orbit on the `i`th section. This last form allows for non-vector state space which can be convenient for 2d problems for example.

- `pb(orbitguess, par)` evaluates the functional G on `orbitguess`
- `pb(orbitguess, par, du)` evaluates the jacobian `dG(orbitguess).du` functional at `orbitguess` on `du`

!!! tip "Tip"
    You can use the function `getPeriod(pb, sol, par)` to get the period of the solution `sol` for the problem with parameters `par`.
"""
@with_kw struct PoincareShootingProblem{Tf, Tsection <: HyperplaneSections} <: AbstractShootingProblem
	M::Int64 = 0								# number of Poincaré sections
	flow::Tf = Flow()							# should be a Flow{TF, Tf, Td}
	section::Tsection = HyperplaneSections(M)	# Poincaré sections
	δ::Float64 = 0e-8							# Numerical value used for the Matrix-Free Jacobian by finite differences. If set to 0, analytical jacobian is used
	isparallel::Bool = false					# whether we use DE in Ensemble mode for multiple shooting
end

@inline isParallel(psh::PoincareShootingProblem) = psh.isparallel

function PoincareShootingProblem(F, p,
			prob::ODEProblem, alg,
			hyp::HyperplaneSections;
			δ = 1e-8, interp_points = 50, isparallel = false, kwargs...)
	pSection(out, u, t, integrator) = (hyp(out, u); out .= out .* (integrator.iter > 1))
	affect!(integrator, idx) = terminate!(integrator)
	# we put nothing option to have an upcrossing
	cb = VectorContinuousCallback(pSection, affect!, hyp.M; interp_points = interp_points, affect_neg! = nothing)
	# change the ODEProblem -> EnsembleProblem for the parallel case
	_M = hyp.M
	isparallel = _M==1 ? false : isparallel
	_pb = isparallel ? EnsembleProblem(prob) : prob
	return PoincareShootingProblem(flow = Flow(F, p, _pb, alg; callback = cb, kwargs...), M = hyp.M, section = hyp, δ = δ, isparallel = isparallel)
end

# this is the "simplest" constructor to use in automatic branching from Hopf
# this is a Hack to pass the arguments to construct a Flow. Indeed, we need to provide the
# appropriate callback for Poincare Shooting to work
PoincareShootingProblem(M::Int, par, prob::ODEProblem, alg; isparallel = false, kwargs...) = PoincareShootingProblem(M = M, flow = (par = par, prob = prob, alg = alg), isparallel = isparallel)

PoincareShootingProblem(M::Int, par, prob1::ODEProblem, alg1, prob2::ODEProblem, alg2; isparallel = false, kwargs...) = PoincareShootingProblem(M = M, flow = (par = par, prob1 = prob1, alg1 = alg1, prob2 = prob2, alg2 = alg2), isparallel = isparallel)

function PoincareShootingProblem(F, p,
			prob::ODEProblem, alg,
			normals::AbstractVector, centers::AbstractVector;
			δ = 1e-8, interp_points = 50, isparallel = false, kwargs...)
	return PoincareShootingProblem(F, p,
					prob, alg,
					HyperplaneSections(normals, centers);
					δ = δ, interp_points = interp_points, isparallel = isparallel, kwargs...)
end

function PoincareShootingProblem(F, p,
				prob1::ODEProblem, alg1,
				prob2::ODEProblem, alg2,
				hyp::HyperplaneSections;
				δ = 1e-8, interp_points = 50, isparallel = false, kwargs...)
	pSection(out, u, t, integrator) = (hyp(out, u); out .= out .* (integrator.iter > 1))
	affect!(integrator, idx) = terminate!(integrator)
	# we put nothing option to have an upcrossing
	cb = VectorContinuousCallback(pSection, affect!, hyp.M; interp_points = interp_points, affect_neg! = nothing)
	# change the ODEProblem -> EnsembleProblem for the parallel case
	_M = hyp.M
	isparallel = _M==1 ? false : isparallel
	_pb1 = isparallel ? EnsembleProblem(prob1) : prob1
	_pb2 = isparallel ? EnsembleProblem(prob2) : prob2
	return PoincareShootingProblem(flow = Flow(F, p, _pb1, alg1, _pb2, alg2; callback = cb, kwargs...), M = hyp.M, section = hyp, δ = δ, isparallel = isparallel)
end

function PoincareShootingProblem(F, p,
				prob1::ODEProblem, alg1,
				prob2::ODEProblem, alg2,
				normals::AbstractVector, centers::AbstractVector;
				δ = 1e-8, interp_points = 50, isparallel = false, kwargs...)
	return PoincareShootingProblem(F, p,
					prob1, alg2, prob2, alg2,
					HyperplaneSections(normals, centers);
					δ = δ, interp_points = interp_points, isparallel = isparallel, kwargs...)
end

"""
	update!(pb::PoincareShootingProblem, centers_bar; _norm = norm)

This function updates the normals and centers of the hyperplanes defining the Poincaré sections.
"""
function update!(pb::PoincareShootingProblem, centers_bar; _norm = norm)
	centers = [E(pb.section, centers_bar[ii], ii) for ii in eachindex(centers_bar)]
	normals = [pb.flow.F(c) for c in centers]
	for ii in eachindex(normals)
		normals[ii] ./= _norm(normals[ii])
	end
	update!(pb.section, normals, centers)
end

function getPeriod(psh::PoincareShootingProblem, x_bar, par)

	M = getM(psh)
	Nm1 = div(length(x_bar), M)

	# reshape the period orbit guess
	x_barc = reshape(x_bar, Nm1, M)

	# variable to hold the computed result
	xc = similar(x_bar, Nm1 + 1, M)
	outc = similar(xc)

	Th = eltype(x_barc)
	period = Th(0)

	# we extend the state space to be able to call the flow, so we fill xc
	if ~isParallel(psh)
		for ii in 1:M
			E!(psh.section, view(xc, :, ii), view(x_barc, :, ii), ii)
			# We need the callback to be active here!!!
			period += @views psh.flow(Val(:TimeSol), xc[:, ii], par, Inf64).t
		end
	else
		for ii in 1:M
			E!(psh.section, view(xc, :, ii), view(x_barc, :, ii), ii)
		end
		solOde =  psh.flow(Val(:TimeSol), xc, par, repeat([Inf64],M))
		for ii in 1:M
			period += solOde[ii].t
		end
	end
	return period
end

function _getMax(psh::PoincareShootingProblem, x_bar::AbstractVector, par; ratio = 1)
	# this function extracts the amplitude of the cycle
	M = getM(psh)
	Nm1 = div(length(x_bar), M)

	# reshape the period orbit guess
	x_barc = reshape(x_bar, Nm1, M)
	xc = similar(x_bar, Nm1 + 1, M)

	Th = eltype(x_bar)
	mx = Th(0)

	if ~isParallel(psh)
		for ii in 1:M
			E!(psh.section, view(xc, :, ii), view(x_barc, :, ii), ii)
			# We need the callback to be active here!!!
			sol = @views psh.flow(Val(:Full), xc[:, ii], par, Inf64)
			mx = max(mx, maximum(sol[1:div(Nm1, ratio), :]))
		end
	else
		for ii in 1:M
			E!(psh.section, view(xc, :, ii), view(x_barc, :, ii), ii)
		end
		solOde =  psh.flow(Val(:Full), xc, par, repeat([Inf64], M) )
		for ii in 1:M
			mx = max(mx, maximum(solOde[ii].u[1:div(Nm1, ratio), :]))
		end
	end

	return mx
end

# Poincaré (multiple) shooting with hyperplanes parametrization
function (psh::PoincareShootingProblem)(x_bar::AbstractVector, par; verbose = false)
	M = getM(psh)
	Nm1 = div(length(x_bar), M)

	# reshape the period orbit guess
	x_barc = reshape(x_bar, Nm1, M)

	# TODO the following declaration of xc allocates. It would be better to make it inplace
	xc = similar(x_bar, Nm1 + 1, M)

	# variable to hold the result of the computations
	outc = similar(xc)

	# we extend the state space to be able to call the flow, so we fill xc
	#TODO create the projections on the fly
	for ii in 1:M
		E!(psh.section, view(xc, :, ii), view(x_barc, :, ii), ii)
	end

	if ~isParallel(psh)
		for ii in 1:M
			im1 = (ii == 1 ? M : ii - 1)
			# We need the callback to be active here!!!
			@views outc[:, ii] .= xc[:, ii] .- psh.flow(xc[:, im1], par, Inf64)
		end
	else
		solOde = psh.flow(xc, par, repeat([Inf64],M))
		for ii in 1:M
			im1 = (ii == 1 ? M : ii - 1)
			# We need the callback to be active here!!!
			@views outc[:, ii] .= xc[:, ii] .- solOde[im1][2]
		end
	end

	# build the array to be returned
	out_bar = similar(x_bar)
	out_barc = reshape(out_bar, Nm1, M)
	for i in 1:M
		R!(psh.section, view(out_barc, :, i), view(outc, :, i), i)
	end
	return out_bar
end

function diffPoincareMap(psh::PoincareShootingProblem, x, par, dx, ii::Int)
	normal = psh.section.normals[ii]
	abs(dot(normal, dx)) > 1e-12 && @warn "Vector does not belong to hyperplane!  dot(normal, dx) = $(abs(dot(normal, dx))) and $(dot(dx, dx))"
	# compute the Poincare map from x
	tΣ, solΣ = psh.flow(Val(:TimeSol), x, par, Inf64)
	z = psh.flow.F(solΣ, par)
	# solution of the variational equation at time tΣ
	# We need the callback to be INACTIVE here!!!
	y = psh.flow(x, par, dx, tΣ; callback = nothing).du
	out = y .- (dot(normal, y) / dot(normal, z)) .* z
end

# jacobian of the shooting functional
function (psh::PoincareShootingProblem)(x_bar::AbstractVector, par, dx_bar::AbstractVector)
	δ = psh.δ
	if δ > 0
		# mostly for debugging purposes
		return (psh(x_bar .+  δ .* dx_bar, par) .- psh(x_bar, par)) ./ δ
	end

	# otherwise analytical Jacobian
	M = getM(psh)
	Nm1 = div(length(x_bar), M)

	# reshape the period orbit guess
	x_barc  = reshape( x_bar, Nm1, M)
	dx_barc = reshape(dx_bar, Nm1, M)

	# variable to hold the computed result
	xc  = similar( x_bar, Nm1 + 1, M)
	dxc = similar(dx_bar, Nm1 + 1, M)
	outc = similar(xc)

	# we extend the state space to be able to call the flow, so we fill xc
	for ii in 1:M
		 E!(psh.section,  view(xc, :, ii),  view(x_barc, :, ii), ii)
		dE!(psh.section, view(dxc, :, ii), view(dx_barc, :, ii), ii)
	end

	if ~isParallel(psh)
		for ii in 1:M
			im1 = (ii == 1 ? M : ii - 1)
			@views outc[:, ii] .= dxc[:, ii] .- diffPoincareMap(psh, xc[:, im1], par, dxc[:, im1], im1)
		end
	else
		@assert 1==0 "Analytical Jacobian for parallel Poincare Shooting not implemented yet. Please use the option δ > 0."
	end

	# build the array to be returned
	out_bar = similar(x_bar)
	out_barc = reshape(out_bar, Nm1, M)
	for ii in 1:M
		dR!(psh.section, view(out_barc, :, ii), view(outc, :, ii), ii)
	end
	return out_bar
end

####################################################################################################
# functions needed Branch switching from Hopf bifurcation point

function update(prob::PoincareShootingProblem, F, dF, hopfpt, ζr, M, centers, period)
	# make the section
	normals = [F(u, hopfpt.params) for u in centers]
	for n in normals; n ./= norm(n); end

	@assert ~(prob.flow isa Flow) "Somehow, this method was not called as it should. prob.flow should be a Named Tuple, prob should be constructed with the simple constructor, not yielding a Flow for its flow field."

	# update the problem
	if length(prob.flow) == 3
	probPSh = PoincareShootingProblem(F, prob.flow.par, prob.flow.prob, prob.flow.alg, normals, centers)
	else
		probPSh = PoincareShootingProblem(F, prob.flow.par, prob.flow.prob1, prob.flow.alg1, prob.flow.prob2, prob.flow.alg2, normals, centers)
	end

	# create initial guess. We have to pass it through the projection R
	hyper = probPSh.section
	orbitguess_bar = zeros(length(centers[1])-1, M)
	for ii=1:length(normals)
		orbitguess_bar[:, ii] .= R(hyper, centers[ii], ii)
	end

	return probPSh, vec(orbitguess_bar)
end