text:
# test_request_format.jl
#
# Author: Tom Price
# Date: Dec 2018
#
# This Julia script tests the JSON data structures
# that are supplied in the body of GET requests
#
# *     based on documentation of REST API from juliaapi_new.txt
#
# *     Julia listens on http://localhost:8081 via HttpServer.jl and dispatch.jl
#       and sends responses to the Rails app on http://localhost:3000
#
# *     the relevant controller in the Rails app is
#       chaipcr/web/app/controllers/experiments_controller.rb
#
# *     currently calls are POST-ed to
#       http://127.0.0.1:8081/experiments/#{experiment.id}/standard_curve
#       http://127.0.0.1:8081/experiments/#{experiment.id}/amplification
#       http://127.0.0.1:8081/experiments/#{experiment.id}/meltcurve
#       http://127.0.0.1:8081/experiments/#{experiment.id}/analyze

import JSON, DataStructures.OrderedDict
import QpcrAnalysis.Action_DICT

# ********************************************************************************
#
# call: experiments/:experiment_id/standard_curve
#
# ********************************************************************************

function standard_curve_request_test()
    request=JSON.parse("""[
        {"well": [
            {"target": 1, "cq": 17.9, "quantity": {"m": 31.2, "b": -3.52}},
            {"target": 2, "cq": 18.1, "quantity": {"m": 27.8, "b": -3.61}}
        ]},
        {"well": [
            {"target": 1, "cq": 17.0, "quantity": {"m": 30.1, "b": -3.47}},
            {"target": 2, "cq": 17.8, "quantity": {"m": 27.5, "b": -3.58}}
        ]},
        {"well": [
            {"target": 1, "cq": 16.2, "quantity": {"m": 29.4, "b": -3.51}},
            {"target": 2, "cq": 17.9, "quantity": {"m": 28.2, "b": -3.55}}
        ]}
    ]"""; dicttype=OrderedDict)
    verify_request(
        Action_DICT["standard_curve"](),
        request
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/amplification
#
# ********************************************************************************

# single channel

function singlechannel_amplification_request_test()
    request=JSON.parse("""{
        "experiment_id": 99,
        "step_id": 99,
        "min_reliable_cyc": 5,
        "calibration_info": {
            "water": {
                "fluorescence_value": [
                    [0.01, 0.02,    0.15, 0.16],
                    null
                ]
            },
            "channel_1": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    null
                ]
            }
        },
        "baseline_cyc_bounds": [5, 10],
        "baseline_method": "sigmoid",
        "cq_method": "Cy0",
        "min_fluomax": 4356,
        "min_D1max": 472,
        "min_D2max": 41,
        "raw_data": {
            "fluorescence_value": [],
            "well_num": [],
            "cycle_num": [],
            "channel": []
        }
    }"""; dicttype=OrderedDict)
    verify_request(
        Action_DICT["amplification"](),
        request
    )
end

# dual channel

function dualchannel_amplification_request_test()
    request=JSON.parse("""{
        "experiment_id": 99,
        "ramp_id": 99,
        "min_reliable_cyc": 5,
        "calibration_info": {
            "water": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    [2.01, 2.02,    2.15, 2.16]
                ]
            },
            "channel_1": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    [2.01, 2.02,    2.15, 2.16]
                ]
            },
            "channel_2": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    [2.01, 2.02,    2.15, 2.16]
                ]
            }
        },
        "baseline_cyc_bounds": [],
        "baseline_method": "sigmoid",
        "cq_method": "Cy0",
        "min_fluomax": 4356,
        "min_D1max": 472,
        "min_D2max": 41,
        "raw_data": {
            "fluorescence_value": [],
            "well_num": [],
            "cycle_num": [],
            "channel": []
        }
    }"""; dicttype=OrderedDict)
    verify_request(
        Action_DICT["amplification"](),
        request
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/meltcurve
#
#
# ********************************************************************************

# single channel

function singlechannel_meltcurve_request_test()
    request=JSON.parse("""{
        "experiment_id": 2,
        "stage_id": 5,
        "calibration_info": {
            "water": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    null
                ]
            },
            "channel_1": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    null
                ]
            },
            "channel_2": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    null
                ]
            }
        },
        "channel_nums": [1],
        "qt_prob": 0.64,
        "max_normd_qtv": 0.8,
        "top_N": 4,
        "raw_data": {
            "fluorescence_value": [],
            "temperature": [],
            "well_num": [],
            "channel": []
        }
    }"""; dicttype=OrderedDict)
    verify_request(
        Action_DICT["meltcurve"](),
        request
    )
end

# dual channel

function dualchannel_meltcurve_request_test()
    request=JSON.parse("""{
        "experiment_id": 2,
        "stage_id": 5,
        "calibration_info": {
            "water": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    [2.01, 2.02,    2.15, 2.16]
                ]
            },
            "channel_1": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    [2.01, 2.02,    2.15, 2.16]
                ]
            },
            "channel_2": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    [2.01, 2.02,    2.15, 2.16]
                ]
            }
        },
        "channel_nums": [1,2],
        "qt_prob": 0.64,
        "max_normd_qtv": 0.8,
        "top_N": 4,
        "raw_data": {
        	"fluorescence_value": [],
        	"temperature": [],
        	"well_num": [],
        	"channel": []
        }
    }"""; dicttype=OrderedDict)
    verify_request(
        Action_DICT["meltcurve"](),
        request
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/thermal_performance_diagnostic
#
# ********************************************************************************

function thermal_performance_diagnostic_request_test()
    request=JSON.parse("""{
      "lid_temp": [1,2,3],
      "heat_block_zone_1_temp": [1,2,3],
      "heat_block_zone_2_temp": [1,2,3],
      "elapsed_time": [1,2,3]
      }"""; dicttype=OrderedDict)
    verify_request(
        Action_DICT["thermal_performance_diagnostic"](),
        request
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/thermal_consistency
#
# ********************************************************************************

# single channel

function single_channel_thermal_consistency_request_test()
    request=JSON.parse("""{
        "raw_data": {
            "fluorescence_value":  [],
            "temperature":  [],
            "well_num":  [],
            "channel":  []
        },
        "calibration_info": {
            "water": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    null
                ]
            },
            "channel_1": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    null
                ]
            },
            "channel_2": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    null
                ]
            }
        }
    }"""; dicttype=OrderedDict)
    verify_request(
        Action_DICT["thermal_consistency"](),
        request
    )
end

# dual channel

function dual_channel_thermal_consistency_request_test()
    request=JSON.parse("""{
        "raw_data": {
            "fluorescence_value":  [],
            "temperature":  [],
            "well_num":  [],
            "channel":  []
        },
        "calibration_info": {
            "water": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    [2.01, 2.02,    2.15, 2.16]
                ]
            },
            "channel_1": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    [2.01, 2.02,    2.15, 2.16]
                ]
            },
            "channel_2": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    [2.01, 2.02,    2.15, 2.16]
                ]
            }
        }
    }"""; dicttype=OrderedDict)
    verify_request(
        Action_DICT["thermal_consistency"](),
        request
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/optical_cal
#
# ********************************************************************************

# single channel

function singlechannel_optical_cal_request_test()
    request=JSON.parse("""{
        "calibration_info": {
            "baseline": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    null
                ]
            },
            "water": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    null
                ]
            },
            "channel_1": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    null
                ]
            },
            "channel_2": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    null
                ]
            }
        }
    }"""; dicttype=OrderedDict)
    verify_request(
        Action_DICT["optical_cal"](),
        request
    )
end

# dual channel

function dualchannel_optical_cal_request_test()
    request=JSON.parse("""{
        "calibration_info": {
            "baseline": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    [2.01, 2.02,    2.15, 2.16]
                ]
            },
            "water": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    [2.01, 2.02,    2.15, 2.16]
                ]
            },
            "channel_1": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    [2.01, 2.02,    2.15, 2.16]
                ]
            },
            "channel_2": {
                "fluorescence_value": [
                    [1.01, 1.02,    1.15, 1.16],
                    [2.01, 2.02,    2.15, 2.16]
                ]
            }
        }
    }"""; dicttype=OrderedDict)
    verify_request(
        Action_DICT["optical_cal"](),
        request
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/optical_test_single_channel
#
# ********************************************************************************

function singlechannel_optical_request_test()
    request=JSON.parse("""{
        "baseline": {
            "fluorescence_value": [1.01, 1.02,    1.15, 1.16]
        },
        "excitation": {
            "fluorescence_value": [1.01, 1.02,    1.15, 1.16]
        }
    }"""; dicttype=OrderedDict)
    verify_request(
        Action_DICT["optical_test_single_channel"](),
        request
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/optical_test_dual_channel
#
# ********************************************************************************

function dualchannel_optical_request_test()
    request=JSON.parse("""{
        "baseline": {
            "fluorescence_value": [
                [1.01, 1.02,    1.15, 1.16],
                [2.01, 2.02,    2.15, 2.16]
            ]
        },
        "water": {
            "fluorescence_value": [
                [1.01, 1.02,    1.15, 1.16],
                [2.01, 2.02,    2.15, 2.16]
            ]
        },
        "channel_1": {
            "fluorescence_value": [
                [1.01, 1.02,    1.15, 1.16],
                [2.01, 2.02,    2.15, 2.16]
            ]
        },
        "channel_2": {
            "fluorescence_value": [
                [1.01, 1.02,    1.15, 1.16],
                [2.01, 2.02,    2.15, 2.16]
            ]
        }
    }"""; dicttype=OrderedDict)
    verify_request(
        Action_DICT["optical_test_dual_channel"](),
        request
    )
end

# ********************************************************************************
#
#                               R U N   E X A M P L E S
#
# ********************************************************************************

function verify_request_examples()
    examples = [
        :standard_curve_request_test,
        :singlechannel_amplification_request_test,
        :dualchannel_amplification_request_test,
        :singlechannel_meltcurve_request_test,
        :dualchannel_meltcurve_request_test,
        :thermal_performance_diagnostic_request_test,
        :single_channel_thermal_consistency_request_test,
        :dual_channel_thermal_consistency_request_test,
        :singlechannel_optical_cal_request_test,
        :dualchannel_optical_cal_request_test,
        :singlechannel_optical_request_test,
        :dualchannel_optical_request_test
    ]
    FactCheck.clear_results()
    OrderedDict(map(examples) do f
        f =>
            try
                !getfield(Main,f)()
            catch err
                err
            end
    end)
end

# Usage:
# verify_request_examples() # every test should return false

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
struct NaiveBayes
    observed_class_distribution::Dict{Int, Float64}
    attribute_observers::Dict{Int, NumericAttributeObserver}
    classes::Vector{Int}
end

function NaiveBayes()
    return NaiveBayes(Dict(), Dict(), [])
end

function partial_fit!(c::NaiveBayes, X::Vector{Float64}, y::Int64, weight::Float64=1.0)
    if haskey(c.observed_class_distribution, y)
        c.observed_class_distribution[y] += weight
    else
        c.observed_class_distribution[y] = weight
    end
    for i in 1:length(X)
        if !haskey(c.attribute_observers, i)
            c.attribute_observers[i] = NumericAttributeObserver()
        end
        update(c.attribute_observers[i], X[i], y)
    end

end

function argmax_votes(votes::Dict{Int64, Float64})
    maxkey = 0
    maxvalue = 0
    for (key, value) in votes
        if value >= maxvalue
            maxkey = key
            maxvalue = value
        end
    end
    return maxkey
end

function predict_proba(c::NaiveBayes, X::Vector{Float64})
    sum_seen_weights = sum(collect(values(c.observed_class_distribution)))
    votes = Dict(i => c.observed_class_distribution[i]/sum_seen_weights for i in keys(c.observed_class_distribution))
    for feature_idx in keys(c.attribute_observers)
        for class_idx in keys(votes)
            class_likelihood = pdf(c.attribute_observers[feature_idx].distribution_per_class[class_idx], X[feature_idx])
            votes[class_idx] *= class_likelihood
        end
    end
    return votes
end

function predict(c::NaiveBayes, X::Vector{Float64})
    votes = predict_proba(c, X)
    return argmax_votes(votes)
end

ext: jl


text:
# Autogenerated wrapper script for gsasl_jll for x86_64-linux-gnu
export gsasl, libgsasl

JLLWrappers.@generate_wrapper_header("gsasl")
JLLWrappers.@declare_executable_product(gsasl)
JLLWrappers.@declare_library_product(libgsasl, "libgsasl.so.18")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        gsasl,
        "bin/gsasl",
    )

    JLLWrappers.@init_library_product(
        libgsasl,
        "lib/libgsasl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    `time_sec(run)` calculates the time in sec.
"""
function time_sec(run::Run)
    cumsum(diff(run.date) / Base.Dates.Millisecond(1000))
end

"""
    `dist_m(run, smoothing = true, λ = 500.0)` estimation of the distance in meters.
"""
function dist_m(run::Run;
    smoothing::Bool = true, λ::AbstractFloat = 500.0)
    lon = run.lon * pi / 180
    lat = run.lat * pi / 180
    m1 = Vector{Float64}(length(lat) - 1)
    for i in 2:length(lat)
        x = (lon[i] - lon[i - 1]) * cos(mean(lat[[i - 1, i]]))
        y = lat[i] - lat[i - 1]
        # authalic radius of 6371007.2 meters
        m1[i - 1] = sqrt(x^2 + y^2) * 6371007.2
    end
    if smoothing
        sec = time_sec(run)
        m1 = predict(fit(SmoothingSpline, sec, m1, λ))
        m1[find(m1 .< 0.0)] .= 0
    end
    m = cumsum(m1)
    m
end

"""
    `speed_kmh(run, smoothing = true, λ = 500.0)` calculates the speed in km/h.
"""
function speed_kmh(run::Run;
    smoothing::Bool = true, λ::AbstractFloat = 500.0)
    m1 = diff(vcat(0.0, dist_m(run, smoothing = smoothing, λ = λ)))
    sec = time_sec(run)
    sec1 = diff(vcat(0.0, sec))
    kmh = (m1 ./ 1000) ./ (sec1 ./ 3600)
    kmh
end

"""
    `pace_minkm(run, smoothing = true, λ = 500.0)` calculates the pace in min/km.
"""
function pace_minkm(run::Run;
    smoothing::Bool = true, λ::AbstractFloat = 500.0)
    m = diff(vcat(0.0, dist_m(run, smoothing = smoothing, λ = λ)))
    sec = diff(vcat(0.0, time_sec(run)))
    minkm = (sec ./ 60) ./ (m ./ 1000)
    minkm
end


text:
using PyCall
pygui(:qt)
using PyPlot
matplotlib[:rcdefaults]()
ion()

using Unwrap
using DSP

x = linspace(0,10,3000)*2π
φ = angle(exp(im*(x-8π).^2))
φφ = @time unwrap1d(φ)
φφφ = @time unwrap(φ)

figure("1d",figsize=(5,4))
clf()

subplot(211)
plot(x/2π,φ/π)
margins(0,0.1)
gca()[:set_xticklabels]([])
ylabel(L"Wrapped [$\pi$]")
subplot(212)
plot(x/2π,φφ/π,label="Unwrap.jl")
plot(x/2π,φφφ/π,label="DSP.jl")
margins(0,0.1)
ylabel(L"Unwrapped [$\pi$]")
legend(framealpha=0.75)
tight_layout()
savefig("1d.png",dpi=72,transparent=true)


text:
using LidJul
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using Makie
using IterativeSolvers
using Preconditioners
using Random

function addmeasurements(msm,solvername,ti,ts,res,iterations=[])
    push!(msm,(solver=solvername,init_time=ti,solver_time=ts,residual=res,iterations=iterations))
end

function resnorm(splu,pref,b)
    ax=similar(b)
    LidJul.mul!(ax,splu,pref)
    norm(ax-b)
end

function test_solver!(msm,Solver,solver_args,splu,b)
    p=zeros(size(b))
    init_time=@elapsed solver=Solver(solver_args...)
    init_time=@elapsed solver=Solver(solver_args...)
    Random.seed!(1234)
    prand=rand(length(p))
    iterations=begin copyto!(p,prand) ; solve!(p,b,solver)  end
    solve_time=@elapsed  begin copyto!(p,prand) ; solve!(p,b,solver) end
    sname=string(Solver)
    length(solver_args)>1 && (sname=sname*"_"*string(solver_args[2]))

    addmeasurements(msm,sname,init_time,solve_time,resnorm(splu,p,b),iterations)
    p
end

function precoconstructor(A,preco)
    if preco==CholeskyPreconditioner
        return preco(A,4)
    elseif preco===nothing
        return IterativeSolvers.Identity()
    else
        return preco(A)
    end
end

include("measureplots.jl")

function testpoisson_GMG(n,bc)
    msm=Vector{NamedTuple{(:solver,:init_time,:solver_time,:residual,:iterations),Tuple{String,Float64,Float64,Float64,Any}}}()

    p=ones(n,n)
    pref=ones(n,n)
    b=zeros(n,n)
    #LP: Create a Poisson's RHS compatible with a possible full Neumann boundary conditions
    #See bellow
    for i=1:n
        for j=1:n
            b[i,j]=sin(3π*(i-1)/(n-1))*sin(6π*(j-1)/(n-1))
        end
    end

    L=1.0
    dx=L/n
    dy=L/n
    @show dx,dy

    l2D=Laplacian2D(n,n,L,L,bc...)
    if bc==(neumann,neumann,neumann,neumann)
        #LP: for all neumann conditions the laplacian equation has no
        #solution if the rhs b does not belong to the range of A.
        #This impose the integral of b to be zero on the domain and
        # the sum of b on the mesh to be also zero...
        #See Mc Cormick and other for more details.
        @show sum(b)
        @assert(sum(b)<1.e-8)

        #Note that if this compatibility condition is ensured, and
        # x is a solution, then x+Cte is also a solution.
    end

    @time s2D=sparse_corr(l2D)
    splu,tlu=PoissonSparseLU(s2D)
    ts=@elapsed solve!(pref,b,splu)

    @time s2D_nocorr=SparseArrays.sparse(l2D)
    splu_nocorr,tlu_nocorr=PoissonSparseLU(s2D_nocorr)
    test_solver!(msm,PoissonGMG,(l2D,GSSmoother),splu_nocorr,b)

    test_solver!(msm,PoissonGMG_new,(l2D,GSSmoother),splu_nocorr,b)

    sname=string(n)*tostring(bc)

    @save "msm"*sname*".jld2" msm

    # measureplot(msm,sname)

    return p,pref,p,splu,b,msm

end

function tostring(bcs)
     r=string()
     for bc in bcs
         bc==dirichlet && (r*="D")
         bc==neumann && (r*="N")
    end
    r
end

function dump_measurments(msm,sname)
    header=keys(msm[1])[1:end]
    nrows=length(msm)
    ncols=length(header)
    table=Matrix{Any}(undef,nrows,ncols)
    iterations=Vector{Union{Nothing,Vector{Float64}}}(undef,nrows)
    for i=1:nrows
        for j=1:ncols
            table[i,j]=msm[i][header[j]]
        end
        table[i]=replace(table[i],"Poisson"=>"")
        table[i]=replace(table[i],"Preconditioner"=>"")
        table[i]=replace(table[i],"SparseCG"=>"PCG")
    end
     #Remove some measurements :
     removed=["GMG_CGSmoother"]
     table=table[filter(i->!(table[i,1] in removed),1:size(table,1)),:]

     #add Total Time and niter
     header=[header[1:3]...,"Total Time",header[4],"niters",header[5]]
     table=hcat(table[:,1:3],table[:,2].+table[:,3],table[:,4],niterations.(residuals(table[:,5])),table[:,5])
     #Sort
     sp=sortperm(table[:, 4])
     table=table[sp,:]
     # iterations=iterations[sp,:]
     nr,nc=size(table)
     println("all methods")
     pretty_table(table[:,1:6],collect(header[1:6]);formatters=ft_printf("%5.3E",2:5),alignment=:l)
end

function go()
    #Choose a power of two
    n=512
    #Choose boundary conditions
    # bc=(neumann,neumann,neumann,neumann)
    bc=(dirichlet,neumann,dirichlet,neumann)
    # bc=(dirichlet,dirichlet,dirichlet,dirichlet)

    sname=string(n)*tostring(bc)
    p,pref,pbj,s,b,msm=testpoisson_GMG(n,bc)
    @save "msm"*sname*".jld2" msm

    dump_measurments(msm,sname)
    @show msm
    msm

   #*** Makie Plots
    x = range(0,1,length=n)
    y = range(0,1,length=n)

    p=pref

    pmax=maximum(p)-minimum(p)

    p./=pmax

    s=Makie.surface(x,y,p)
    xm, ym, zm = minimum(scene_limits(s))
    Makie.contour!(s,x,y,p, levels = 15, linewidth = 2, transformation = (:xy, zm))
    Makie.wireframe!(s,x,y,p, overdraw = true, transparency = true, color = (:black, 0.1))
    display(AbstractPlotting.PlotDisplay(), s)
    resize!(s,(1600,800))
    Makie.save("makie"*string(n)*tostring(bc)*".png", s)
    nothing
end

msm=go()


text:
using Revise

using JDF, CSV, DataFrames, Blosc, Base.GC, BenchmarkTools, Serialization

# use 12 threads
Blosc.set_num_threads(6)
# GC.gc()
@time a = CSV.read(
    "C:/data/Performance_All/Performance_2010Q3.txt",
    delim = '|',
    header = false
);
GC.gc()

@time metadatas = savejdf("c:/data/large.jdf", a);

using Revise, JDF, DataFrames, Serialization
@time a2 = loadjdf("c:/data/large.jdf"; verbose);
GC.gc()

@time type_compress!(a2);

io = iow()
compress_then_write(a2.Column2, io)
close(io)

@time psavejdf("c:/data/largec.jdf", a2);

using Revise, JDF, DataFrames, Serialization
@time a2 = ploadjdf("c:/data/largec.jdf");
GC.gc()

@time type_compress!(a2, verbose = true)

@time psavejdf("c:/data/largec2.jdf", a2);
@time a2 = loadjdf("c:/data/largec.jdf");

using Revise, JDF, DataFrames, Serialization
@time a2 = ploadjdf("c:/data/largec.jdf");
GC.gc()

all(names(a) .== names(a2))
all(skipmissing([all(a2[!,name] .== Array(a[!,name])) for name in names(a2)]))

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
#
# Solving the SPMe and comparison with the PyBaMM solution
#

using PyBaMM
using SparseArrays, LinearAlgebra
using PyCall

# load model
pybamm = pyimport("pybamm")
model = pybamm.lithium_ion.SPMe(name="SPMe")
sim = pybamm.Simulation(model)

prob = get_ode_problem(sim)

using Sundials
sol = solve(prob, CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1), reltol=1e-6, abstol=1e-6, saveat=prob.tspan[2] / 100);

# Calculate voltage in Julia
V = get_variable(sim, sol, "Terminal voltage [V]")
t = get_variable(sim, sol, "Time [s]")

# Solve in python
sol_pybamm = sim.solve(t)
V_pybamm = get(sol_pybamm, "Terminal voltage [V]").data

# Plots
using Plots

plot(sol.t, V)
scatter!(sol.t, V_pybamm)


text:
# Tests for Reverse DDF DEA Models
@testset "ReverseDDFDEAModel" begin

    # ------------------
    # Input oriented
    # ------------------
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6 8]
    Y = [1; 1; 1; 1; 1; 1; 1; 1]

    # Reverse DDF for input :ERG
    russelio = dearussell(X, Y, orient = :Input, rts = :VRS)
    rddfergio = dearddf(X, Y, :ERG, orient = :Input, rts = :VRS)

    @test efficiency(rddfergio) ≈ [0.0; 0.0; 0.0; 5/12; 0.6; 1/6; 0.35; 0.4375] atol = 1e-5
    @test targets(russelio, :X) ≈ targets(rddfergio, :X)
    @test targets(russelio, :Y) ≈ targets(rddfergio, :Y)

    # Print
    show(IOBuffer(), rddfergio)

    # ------------------
    # Output oriented
    # ------------------
    X = [1; 1; 1; 1; 1; 1; 1; 1]
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5]

    # Reverse DDF for output :ERG
    russeloo = dearussell(X, Y, orient = :Output, rts = :VRS)
    rddfergoo = dearddf(X, Y, :ERG, orient = :Output, rts = :VRS)

    @test efficiency(rddfergoo) ≈ [0.0; 0.0; 0.0; 13/28; 4/7; 1/3; 11/35; 37/55] atol = 1e-5
    @test targets(russeloo, :X) ≈ targets(rddfergoo, :X)
    @test targets(russeloo, :Y) ≈ targets(rddfergoo, :Y)

    # Print
    show(IOBuffer(), rddfergoo)

    # ------------------
    # Graph
    # ------------------
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]

    # Reverse DDF for :ERG :VRS
    erg = deaerg(X, Y, rts = :VRS)
    rddferg = dearddf(X, Y, :ERG, orient = :Graph, rts = :VRS)

    @test efficiency(rddferg) ≈ [0.0; 0.0; 0.0; 0.0; 0.6; 11/21; 1/7; 0.8] atol = 1e-5
    @test targets(erg, :X) ≈ targets(rddferg, :X)
    @test targets(erg, :Y) ≈ targets(rddferg, :Y)

    # Reverse DDF for :ERG :VCRS
    ergcrs = deaerg(X, Y, rts = :VRS)
    rddfergcrs = dearddf(X, Y, :ERG, orient = :Graph, rts = :VRS)

    @test targets(ergcrs, :X) ≈ targets(rddfergcrs, :X)
    @test targets(ergcrs, :Y) ≈ targets(rddfergcrs, :Y)

    # Default orient is :Graph
    @test efficiency(dearddf(X, Y, :ERG, rts = :VRS)) == efficiency(rddferg)

    # Print
    show(IOBuffer(), rddferg)

    # Reverse DDF for Modified DDF for :Observed :VRS
    mddf = deamddf(X, Y, Gx = :Observed, Gy = :Observed, rts = :VRS, slack = false)
    rddfmddf = dearddf(X, Y, :MDDF, Gx = :Observed, Gy = :Observed, rts = :VRS)

    @test efficiency(rddfmddf) ≈ [0.0; 0.0; 0.0; 0.0; 7/6; 4/7; 1/7; 2.549936] atol = 1e-5
    @test targets(rddfmddf, :X) ≈ targets(mddf, :X) atol = 1e-5
    @test targets(rddfmddf, :Y) ≈ targets(mddf, :Y) atol = 1e-5

    # Reverse DDF for Modified DDF for :Observed :CRS
    mddfcrs = deamddf(X, Y, Gx = :Observed, Gy = :Observed, rts = :CRS, slack = false)
    rddfmddfcrs = dearddf(X, Y, :MDDF, Gx = :Observed, Gy = :Observed, rts = :CRS)

    @test targets(rddfmddfcrs, :X) ≈ targets(mddfcrs, :X) atol = 1e-5
    @test targets(rddfmddfcrs, :Y) ≈ targets(mddfcrs, :Y) atol = 1e-5

    # Reverse DDF for Modified DDF for :Ones :VRS
    mddfones = deamddf(X, Y, Gx = :Ones, Gy = :Ones, rts = :VRS, slack = false)
    rddfmddfones = dearddf(X, Y, :MDDF, Gx = :Ones, Gy = :Ones, rts = :VRS)

    @test targets(rddfmddfones, :X) ≈ targets(mddfones, :X) atol = 1e-5
    @test targets(rddfmddfones, :Y) ≈ targets(mddfones, :Y) atol = 1e-5

    # Reverse DDF for Modified DDF for :Mean :VRS
    mddfmean = deamddf(X, Y, Gx = :Mean, Gy = :Mean, rts = :VRS, slack = false)
    rddfmddfmean = dearddf(X, Y, :MDDF, Gx = :Mean, Gy = :Mean, rts = :VRS)

    @test targets(rddfmddfmean, :X) ≈ targets(mddfmean, :X) atol = 1e-5
    @test targets(rddfmddfmean, :Y) ≈ targets(mddfmean, :Y) atol = 1e-5

    # Print
    show(IOBuffer(), rddfmddf)

    ## Test if one-by-one DEA using evaluation and reference sets match initial results
    rddf_vrs_ref_eff = zeros(size(X, 1))

    Xref = X[:,:]
    Yref = Y[:,:]

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]

        rddf_vrs_ref_eff[i] = efficiency(dearddf(Xeval, Yeval, :ERG; rts = :VRS, Xref = Xref, Yref = Yref))[1]
    end

    @test rddf_vrs_ref_eff ≈ efficiency(rddferg)

    # ------------------
    # Test slacks with different data
    # ------------------
    X = [1 1; 2 2; 1.5 1]
    Y = [1 1; 2 2; 1 0.5]

    mddfslack = deamddf(X, Y, Gx = :Observed, Gy = :Observed, rts = :VRS)
    rddfslack = dearddf(X, Y, :MDDF, Gx = :Observed, Gy = :Observed, rts = :VRS, slack = true)
    @test slacks(rddfslack, :X) ≈ [0 0; 0 0; 0.5 0] atol = 1e-5
    @test slacks(rddfslack, :Y) ≈ [0 0; 0 0; 0 0.5] atol = 1e-5

    # Print
    show(IOBuffer(), rddfslack)

    # ------------------
    # Test errors
    # ------------------
    @test_throws DimensionMismatch dearddf([1; 2 ; 3], [4 ; 5], :ERG) #  Different number of observations
    @test_throws DimensionMismatch dearddf([1; 2], [4 ; 5], :ERG, Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws DimensionMismatch dearddf([1 1; 2 2], [4 4; 5 5], :ERG, Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws DimensionMismatch dearddf([1 1; 2 2], [4 4; 5 5], :ERG, Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ArgumentError dearddf([1; 2; 3], [4; 5; 6], :Error) # Invalid efficiency measure
    @test_throws ArgumentError dearddf([1; 2; 3], [4; 5; 6], :ERG, rts = :Error) # Invalid returns to scale
    @test_throws ArgumentError dearddf([1; 2; 3], [4; 5; 6], :ERG, orient = :Error) # Invalid orientation
    @test_throws ArgumentError dearddf([1; 2; 3], [4; 5; 6], :ERG, Gx = :Ones) # No direction in :ERG
    @test_throws ArgumentError dearddf([1; 2; 3], [4; 5; 6], :MDDF, orient = :Input, Gx = :Ones, Gy = :Ones) # Invalid direction in :MDDF
    @test_throws DimensionMismatch dearddf([1 1; 2 2; 3 3], [4; 5; 6], :MDDF, Gx = [1 1 1; 2 2 2; 3 3 3], Gy = [4; 5; 6]) # Different size of inputs direction
    @test_throws DimensionMismatch dearddf([1; 2; 3], [4 4; 5 5; 6 6], :MDDF, Gx = [1; 2; 3], Gy = [4 4 4; 5 5 5; 6 6 6]) # Different size of inputs direction
    @test_throws ArgumentError dearddf([1; 2; 3], [1; 2; 3], :MDDF, Gx = :Error, Gy = :Ones) # Invalid inuts direction
    @test_throws ArgumentError dearddf([1; 2; 3], [1; 2; 3], :MDDF, Gx = :Ones, Gy = :Error) # Invalid outputs direction

end


text:
@coeff function K_b(reactor::Reactor)
  nu_term = (1+nu_n)*(1+nu_T) * (nu_n + 0.054*nu_T)
  0.6099 * epsilon^(5/2) * kappa^1.27 * nu_term * C_B
end


text:
function test_complex_1()

    Random.seed!(0)
    f = x -> [exp(-30 / x), exp(-30 / x)]
    xtrue = 100
    y = f(xtrue)
    T = 1
    xDists = Uniform(10, 500)
    noiseDist = Normal(0, 0.01)
    signalModels = f
    λ = 2.0^-1.5
    kernel = GaussianKernel(λ * y)
    ρ = 2.0^-20
    xhat = perk(reshape(y, :, 1), T, xDists, noiseDist, signalModels, kernel, ρ)
    error_rel_real = abs(xhat[] - xtrue) / xtrue

    Random.seed!(0)
    f = x -> complex(exp(-30 / x), exp(-30 / x))
    xtrue = 100
    y = f(xtrue)
    T = 1
    xDists = Uniform(10, 500)
    noiseDist = Normal(0, 0.01)
    signalModels = f
    λ = 2.0^-1.5
    kernel = GaussianKernel([λ * mean(y)])
    ρ = 2.0^-20
    xhat = perk(y, T, xDists, noiseDist, signalModels, kernel, ρ)
    error_rel_complex = abs(xhat[] - xtrue) / xtrue

    return error_rel_complex ≈ error_rel_real

end

function test_complex_2()

    Random.seed!(0)
    f = x -> complex(exp(-30 / x), exp(-30 / x))
    xtrue = 100
    y = f(xtrue)
    T = 200
    xDists = Uniform(10, 500)
    noiseDist = Normal(0, 0.01)
    signalModels = f
    λ = 2.0^-1.5
    kernel = GaussianKernel([λ * mean(y)])
    ρ = 2.0^-20
    xhat = perk(y, T, xDists, noiseDist, signalModels, kernel, ρ)

    error_rel = abs(xhat[] - xtrue) / xtrue
    return isapprox(error_rel, 0.05211190543523628, atol = 1e-6)

end

function test_complex_3()

    Random.seed!(0)
    f = (x, ν) -> [exp(-ν / x), exp(-ν / x)]
    xtrue = 100
    ν = 30
    N = 1
    y = reshape(f(xtrue, ν), :, N)
    T = 1
    xDists = [Uniform(10, 500)]
    νDists = Uniform(30, 30.000000001)
    noiseDist = Normal(0, 0.01)
    signalModels = [f]
    λ = 2.0^-1.5
    kernel = GaussianKernel(λ * vec(y), [λ * mean(ν)])
    ρ = 2.0^-20
    xhat = perk(y, ν, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)
    error_rel_real = abs(xhat[] - xtrue) / xtrue

    Random.seed!(0)
    f = (x, ν) -> complex(exp(-ν / x), exp(-ν / x))
    xtrue = 100
    ν = 30
    N = 1
    y = fill(f(xtrue, ν), 1, N)
    T = 1
    xDists = [Uniform(10, 500)]
    νDists = Uniform(30, 30.000000001)
    noiseDist = Normal(0, 0.01)
    signalModels = [f]
    λ = 2.0^-1.5
    kernel = GaussianKernel([λ * mean(y)], [λ * mean(ν)])
    ρ = 2.0^-20
    xhat = perk(y, ν, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)
    error_rel_complex = abs(xhat[] - xtrue) / xtrue

    return error_rel_complex ≈ error_rel_real

end

function test_complex_4()

    Random.seed!(0)
    f = (x, ν) -> complex(exp(-ν / x), exp(-ν / x))
    xtrue = 100
    ν = 30
    N = 1
    y = fill(f(xtrue, ν), 1, N)
    T = 200
    xDists = [Uniform(10, 500)]
    νDists = Uniform(30, 30.000000001)
    noiseDist = Normal(0, 0.01)
    signalModels = [f]
    λ = 2.0^-1.5
    kernel = GaussianKernel([λ * mean(y)], [λ * mean(ν)])
    ρ = 2.0^-20
    xhat = perk(y, ν, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)

    error_rel = abs(xhat[] - xtrue) / xtrue
    return isapprox(error_rel, 0.03526687067013938, atol = 1e-6)

end

function test_complex_5()

    Random.seed!(0)
    f = x -> complex(exp(-30 / x) + 1, exp(-30 / x) + 1)
    f_abs = x -> abs(f(x))
    xtrue = 100
    y = f(xtrue)
    y_abs = f_abs(xtrue)
    T = 200
    xDists = Uniform(10, 500)
    noiseDist = Normal(0, 0.01)
    λ = 2.0^-1.5
    kernel = GaussianKernel([λ * mean(y)])
    kernel_abs = GaussianKernel(λ * mean(y_abs))
    ρ = 2.0^-20
    xhat = perk(y, T, xDists, noiseDist, f, kernel, ρ)
    xhat_abs = perk(y_abs, T, xDists, noiseDist, f_abs, kernel_abs, ρ)

    error_rel = abs(xhat[] - xtrue) / xtrue
    error_rel_abs = abs(xhat_abs[] - xtrue) / xtrue
    return error_rel < error_rel_abs

end

function test_complex_6()

    Random.seed!(0)
    f = x -> [exp(-30 / x), exp(-30 / x)]
    xtrue = 100
    y = f(xtrue)
    T = 1
    xDists = Uniform(10, 500)
    noiseDist = Normal(0, 0.01)
    signalModels = f
    λ = 2.0^-1.5
    H = 2
    kernel = GaussianRFF(H, λ * y)
    ρ = 2.0^-20
    xhat = perk(reshape(y, :, 1), T, xDists, noiseDist, signalModels, kernel, ρ)
    error_rel_real = abs(xhat[] - xtrue) / xtrue

    Random.seed!(0)
    f = x -> complex(exp(-30 / x), exp(-30 / x))
    xtrue = 100
    y = f(xtrue)
    T = 1
    xDists = Uniform(10, 500)
    noiseDist = Normal(0, 0.01)
    signalModels = f
    λ = 2.0^-1.5
    H = 2
    kernel = GaussianRFF(H, [λ * mean(y)])
    ρ = 2.0^-20
    xhat = perk(y, T, xDists, noiseDist, signalModels, kernel, ρ)
    error_rel_complex = abs(xhat[] - xtrue) / xtrue

    return error_rel_complex ≈ error_rel_real

end

function test_complex_7()

    Random.seed!(0)
    f = x -> complex(exp(-30 / x), exp(-30 / x))
    xtrue = 100
    y = f(xtrue)
    T = 200
    xDists = Uniform(10, 500)
    noiseDist = Normal(0, 0.01)
    signalModels = f
    λ = 2.0^-1.5
    H = 40
    kernel = GaussianRFF(H, [λ * mean(y)])
    ρ = 2.0^-20
    xhat = perk(y, T, xDists, noiseDist, signalModels, kernel, ρ)

    error_rel = abs(xhat[] - xtrue) / xtrue
    return isapprox(error_rel, 0.01699497197550329, atol = 1e-2)

end

function test_complex_8()

    Random.seed!(0)
    f = (x, ν) -> [exp(-ν / x), exp(-ν / x)]
    xtrue = 100
    ν = 30
    N = 1
    y = reshape(f(xtrue, ν), :, N)
    T = 1
    xDists = [Uniform(10, 500)]
    νDists = Uniform(30, 30.000000001)
    noiseDist = Normal(0, 0.01)
    signalModels = [f]
    λ = 2.0^-1.5
    H = 2
    kernel = GaussianRFF(H, λ * vec(y), [λ * mean(ν)])
    ρ = 2.0^-20
    xhat = perk(y, ν, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)
    error_rel_real = abs(xhat[] - xtrue) / xtrue

    Random.seed!(0)
    f = (x, ν) -> complex(exp(-ν / x), exp(-ν / x))
    xtrue = 100
    ν = 30
    N = 1
    y = fill(f(xtrue, ν), 1, N)
    T = 1
    xDists = [Uniform(10, 500)]
    νDists = Uniform(30, 30.000000001)
    noiseDist = Normal(0, 0.01)
    signalModels = [f]
    λ = 2.0^-1.5
    H = 2
    kernel = GaussianRFF(H, [λ * mean(y)], [λ * mean(ν)])
    ρ = 2.0^-20
    xhat = perk(y, ν, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)
    error_rel_complex = abs(xhat[] - xtrue) / xtrue

    return error_rel_complex ≈ error_rel_real

end

function test_complex_9()

    Random.seed!(0)
    f = (x, ν) -> complex(exp(-ν / x), exp(-ν / x))
    xtrue = 100
    ν = 30
    N = 1
    y = fill(f(xtrue, ν), 1, N)
    T = 200
    xDists = [Uniform(10, 500)]
    νDists = Uniform(30, 30.000000001)
    noiseDist = Normal(0, 0.01)
    signalModels = [f]
    λ = 2.0^-1.5
    H = 40
    kernel = GaussianRFF(H, [λ * mean(y)], [λ * mean(ν)])
    ρ = 2.0^-20
    xhat = perk(y, ν, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)

    error_rel = abs(xhat[] - xtrue) / xtrue
    return isapprox(error_rel, 0.0346180252991536, atol = 1e-2)

end

function test_complex_10()

    Random.seed!(0)
    f = x -> complex(exp(-30 / x) + 1, exp(-30 / x) + 1)
    f_abs = x -> abs(f(x))
    xtrue = 100
    y = f(xtrue)
    y_abs = f_abs(xtrue)
    T = 200
    xDists = Uniform(10, 500)
    noiseDist = Normal(0, 0.01)
    λ = 2.0^-1.5
    H = 40
    kernel = GaussianRFF(H, [λ * mean(y)])
    kernel_abs = GaussianRFF(H, λ * mean(y_abs))
    ρ = 2.0^-20
    xhat = perk(y, T, xDists, noiseDist, f, kernel, ρ)
    xhat_abs = perk(y_abs, T, xDists, noiseDist, f_abs, kernel_abs, ρ)

    error_rel = abs(xhat[] - xtrue) / xtrue
    error_rel_abs = abs(xhat_abs[] - xtrue) / xtrue
    return error_rel < error_rel_abs

end

@testset "Complex PERK" begin

    @test test_complex_1()
    @test test_complex_2()
    @test test_complex_3()
    @test test_complex_4()
    @test test_complex_5()
    @test test_complex_6()
    @test test_complex_7()
    @test test_complex_8()
    @test test_complex_9()
    @test test_complex_10()

end


text:
"""
**Expected number of outgoing degrees**

    degree_out(N::AbstractEcologicalNetwork)
"""
function degree_out(N::AbstractEcologicalNetwork)
  d_o = vec(sum(N.A, dims=2))
  return Dict(zip(species(N; dims=1), d_o))
end

"""
**Expected number of ingoing degrees**

    degree_in(N::AbstractEcologicalNetwork)
"""
function degree_in(N::AbstractEcologicalNetwork)
  d_i = vec(sum(N.A, dims=1))
  return Dict(zip(species(N; dims=2), d_i))
end

function degree(N::AbstractEcologicalNetwork; dims::Union{Nothing,Integer}=nothing)
  dims == 1 && return degree_out(N)
  dims == 2 && return degree_in(N)
  if dims === nothing
    if typeof(N) <: AbstractBipartiteNetwork
      return merge(degree_out(N), degree_in(N))
    else
      din = degree_in(N)
      dout = degree_out(N)
      for k in keys(din)
        din[k] += dout[k]
      end
      return din
    end
  end
  throw(ArgumentError("dims can only be 1 (out-degre) or 2 (in-degree) or `nothing` (both), you used $(dims)"))
end

function degree(N::QuantitativeNetwork; dims::Union{Nothing,Integer}=nothing)
  degree(convert(BinaryNetwork, N), dims=dims)
end

"""
**Variance in the outgoing degree**

    degree_out_var(N::ProbabilisticNetwork)
"""
function degree_out_var(N::ProbabilisticNetwork)
  d_o_v = mapslices(a_var, N.A; dims=2)
  return Dict(zip(species(N; dims=1), d_o_v))
end

"""
**Variance in the ingoing degree**

    degree_in_var(N::ProbabilisticNetwork)
"""
function degree_in_var(N::ProbabilisticNetwork)
  d_i_v = mapslices(a_var, N.A; dims=1)'
  return Dict(zip(species(N; dims=2), d_i_v))
end

"""
    degree_var(N::ProbabilisticNetwork; dims::Union{Nothing,Integer}=nothing)

Variance in the degree of species in a probabilistic network.
"""
function degree_var(N::ProbabilisticNetwork; dims::Union{Nothing,Integer}=nothing)
  dims == 1 && return degree_out_var(N)
  dims == 2 && return degree_in_var(N)
  if dims === nothing
    if typeof(N) <: AbstractBipartiteNetwork
      return merge(degree_out_var(N), degree_in_var(N))
    else
      din = degree_in_var(N)
      dout = degree_out_var(N)
      for k in keys(din)
        din[k] += dout[k]
      end
      return din
    end
  end
  throw(ArgumentError("dims can only be 1 (out-degre) or 2 (in-degree) or `nothing` (both), you used $(dims)"))
end


text:
"""
    ShapeColl{S<:AbsShapeDesc}

This is the type of a homogeneous collection of finite element shapes.

- `S` = shape descriptor: subtype of AbsShapeDesc.
"""
struct ShapeColl{S<:AbsShapeDesc}
    shapedesc::S # Shape descriptor
    nshapes::Int64 # Number of shapes in the collection
    attributes::Dict{String,AbsAttrib} # dictionary of attributes
    name::String # name of the shape collection
end

"""
    ShapeColl(shapedesc::S, nshapes::Int64) where {S<:AbsShapeDesc}

Set up new shape collection with an empty dictionary of attributes and a default name.
"""
function ShapeColl(shapedesc::S, nshapes::Int64) where {S<:AbsShapeDesc}
    ShapeColl(shapedesc::S, nshapes, Dict{String,AbsAttrib}(), "shapes")
end

"""
    ShapeColl(shapedesc::S, nshapes::Int64, s::String) where {S<:AbsShapeDesc}

Set up new shape collection with an empty dictionary of attributes.
"""
function ShapeColl(shapedesc::S, nshapes::Int64, s::String) where {S<:AbsShapeDesc}
    ShapeColl(shapedesc::S, nshapes, Dict{String,AbsAttrib}(), s)
end

"""
    shapedesc(shapes::ShapeColl)

Retrieve the shape descriptor.
"""
shapedesc(shapes::ShapeColl) = shapes.shapedesc

"""
    nshapes(shapes::ShapeColl)

Number of shapes in the collection.
"""
nshapes(shapes::ShapeColl) = shapes.nshapes

"""
    manifdim(shapes::ShapeColl)

Retrieve the manifold dimension of the collection.
"""
manifdim(shapes::ShapeColl) = manifdim(shapes.shapedesc)

"""
    nvertices(shapes::ShapeColl)

Retrieve the number of vertices per shape.
"""
nvertices(shapes::ShapeColl) = nvertices(shapes.shapedesc)

"""
    facetdesc(shapes::ShapeColl)

Retrieve the shape type of the boundary facet.
"""
facetdesc(shapes::ShapeColl) = shapes.shapedesc.facetdesc

"""
    nfacets(shapes::ShapeColl)

Retrieve the number of boundary facets per shape.
"""
nfacets(shapes::ShapeColl) = nfacets(shapes.shapedesc)

"""
    facetconnectivity(shapes::ShapeColl, j::Int64)

Retrieve connectivity of the `j`-th facet shape of the `i`-th shape from the collection.
"""
function facetconnectivity(shapes::ShapeColl, j::Int64)
    return shapes.shapedesc.facets[j, :]
end

"""
    ridgedesc(shapes::ShapeColl)

Retrieve the shape type of the boundary ridge.
"""
ridgedesc(shapes::ShapeColl) = shapes.shapedesc.ridgedesc

"""
    nridges(shapes::ShapeColl)

Retrieve the number of boundary ridges per shape.
"""
nridges(shapes::ShapeColl) = nridges(shapes.shapedesc)

"""
    ridgeconnectivity(shapes::ShapeColl, j::Int64)

Retrieve connectivity of the `j`-th ridge shape of the `i`-th shape from the collection.
"""
function ridgeconnectivity(shapes::ShapeColl, j::Int64)
    return shapes.shapedesc.ridges[j, :]
end

"""
    attribute(shapes::ShapeColl, s::Symbol)

Retrieve a named attribute.
"""
attribute(shapes::ShapeColl, s::String) = shapes.attributes[s]

function _sense(fc, oc, nshifts) # is the facet used in the positive or in the negative sense?
    if fc == oc
        return +1 # facet used in the positive sense
    end
    for i in 1:(nshifts - 1)
        fc = circshift(fc, 1) # try a circular shift
        if fc == oc
            return +1 # facet used in the positive sense
        end
    end
    return -1 # facet used in the positive sense
end

"""
    Base.summary(sc::S) where {S<:ShapeColl}

Form a brief summary of the shape collection.
"""
function Base.summary(sc::S) where {S<:ShapeColl}
    s = "$(sc.name) = $(nshapes(sc)) x $(shapedesc(sc).name)"
    if !isempty(values(sc.attributes))
        s = s * " {"
        for k in keys(sc.attributes)
            s = s * k * ","
        end
        s = s * "}"
    end

    return s
end


text:
export minnorm

function minnorm(L::Tx, b, k, eta, rank1, iter)where {ETx<:Number, Tx<:AbstractMatrix{ETx}}
    ######Use closed form solution to solve for k-th row of R or L
    bsub = b[:,k];
    ind = findall(x->x!=0, bsub);
    v = minL2NormSPG(L[ind,:],bsub[ind,:],rank1, eta, iter);
    return v
end


text:
#=
    compute the nth root of x
    follows mpfun90.f mpnrt

    iteration converges to a^(-1/n)
    the reciprocal of the final approximation is the nth root

    x[k+1] = x[k] + ((x[k] / n) * (1 - (a * (x[k]^n))

    iterations support lesser precision early, roughly doubling
    the precision with each next stage of the iterative process

=#
function rootn_ddsi_dd(a::Tuple{T,T}, n::Signed) where {T<:AbstractFloat}
    one1 = one(T)
    hi, lo = a
    n_fp = T(n)
    n_inv = inv_fp_dd(n_fp)
    # initial approximation
    x = inv(hi)^(HI(n_inv))

    m = n
    xon = dvi_fpfp_dd(x, n_fp)
    xpn = mul_fpfp_dd(x, x)
    m -= 2
    while m > 0
       xpn = mul_fpdd_dd(x, xpn)
       m -= 1
    end
    axn = mul_dddd_dd(a, xpn)
    axn = sub_fpdd_dd(one1, axn)
    xon = mul_dddd_dd(xon, axn)
    x = add_fpdd_dd(x, xon)

    m = n
    xon = dvi_ddfp_dd(x, n_fp)
    xpn = mul_dddd_dd(x, x)
    m -= 2
    while m > 0
       xpn = mul_dddd_dd(x, xpn)
       m -= 1
    end
    axn = mul_dddd_dd(a, xpn)
    axn = sub_fpdd_dd(one1, axn)
    xon = mul_dddd_dd(xon, axn)
    x = add_dddd_dd(x, xon)

    m = n
    xon = dvi_ddfp_dd(x, n_fp)
    xpn = mul_dddd_dd(x, x)
    m -= 2
    while m > 0
       xpn = mul_dddd_dd(x, xpn)
       m -= 1
    end
    axn = mul_dddd_dd(a, xpn)
    axn = sub_fpdd_dd(one1, axn)
    xon = mul_dddd_dd(xon, axn)
    x = add_dddd_dd(x, xon)

    x = inv_dd_dd(x)
    return x
end

function root5_ddsi_dd(a::Tuple{T,T}, n::Signed) where {T<:AbstractFloat}
    one1 = one(T)
    hi, lo = a
    n_fp = T(n)
    n_inv = inv_fp_dd(n_fp)
    # initial approximation
    x = inv(hi)^(HI(n_inv))

   # m = n
    xon = dvi_fpfp_dd(x, n_fp)
    xpn = mul_fpfp_dd(x, x)
    xpn = mul_fpdd_dd(x, xpn)
    xpn = mul_fpdd_dd(x, xpn)
    xpn = mul_fpdd_dd(x, xpn)
   # m -= 2
   # while m > 0
   #    xpn = mul_fpdd_dd(x, xpn)
   #    m -= 1
   # end
    axn = mul_dddd_dd(a, xpn)
    axn = sub_fpdd_dd(one1, axn)
    xon = mul_dddd_dd(xon, axn)
    x = add_fpdd_dd(x, xon)

    #m = n
    xon = dvi_ddfp_dd(x, n_fp)
    xpn = mul_dddd_dd(x, x)
    xpn = mul_dddd_dd(x, xpn)
    xpn = mul_dddd_dd(x, xpn)
    xpn = mul_dddd_dd(x, xpn)
   #
    #m -= 2
    #while m > 0
    #   xpn = mul_dddd_dd(x, xpn)
    #   m -= 1
    #end
    axn = mul_dddd_dd(a, xpn)
    axn = sub_fpdd_dd(one1, axn)
    xon = mul_dddd_dd(xon, axn)
    x = add_dddd_dd(x, xon)

   # m = n
    xon = dvi_ddfp_dd(x, n_fp)
    xpn = mul_dddd_dd(x, x)
    xpn = mul_dddd_dd(x, xpn)
    xpn = mul_dddd_dd(x, xpn)
    xpn = mul_dddd_dd(x, xpn)

   # m -= 2
   # while m > 0
   #    xpn = mul_dddd_dd(x, xpn)
   #    m -= 1
   # end
    axn = mul_dddd_dd(a, xpn)
    axn = sub_fpdd_dd(one1, axn)
    xon = mul_dddd_dd(xon, axn)
    x = add_dddd_dd(x, xon)

    x = inv_dd_dd(x)
    return x
end


text:
function sosfilters(input::AbstractMatrix, coeffs::AbstractMatrix, state)

    nFilterChannels = size(coeffs, 1)
    nInputChannels, nSamples = size(input)
    nOutputChannels = max(nFilterChannels, nInputChannels)

    # if !(nFilterChannels == nInputChannels == nOutputChannels)
    #     error("Channel numbers don't match: $nFilterChannels filter, $nInputChannels input and $nOutputChannels output channels.")
    # end

    output = similar(input, (nOutputChannels, nSamples))

    if nInputChannels == 1
        error("1 input channel not implemented.")
    elseif nFilterChannels == 1
        a0 = coeffs[1, 1]
        a1 = coeffs[1, 2]
        a2 = coeffs[1, 3]
        b1 = coeffs[1, 4]
        b2 = coeffs[1, 5]

        for n in 1:nSamples
            for i in 1:nOutputChannels
                inp = input[i, n]
                output[i, n] = a0 * inp + state[i, 1]
                state[i, 1] = a1 * inp - b1 * output[i, n] + state[i, 2]
                state[i, 2] = a2 * inp - b2 * output[i, n]
            end
        end
    else
        for n in 1:nSamples
            for i in 1:nOutputChannels

                a0 = coeffs[i, 1]
                a1 = coeffs[i, 2]
                a2 = coeffs[i, 3]
                b1 = coeffs[i, 4]
                b2 = coeffs[i, 5]

                inp = input[i, n]
                output[i, n] = a0 * inp + state[i, 1]
                state[i, 1] = a1 * inp - b1 * output[i, n] + state[i, 2]
                state[i, 2] = a2 * inp - b2 * output[i, n]
            end
        end
    end

    output, state
end


text:
# https://sudonull.com/post/7909-Julia-Scripts-and-parsing-command-line-arguments

using DelimitedFiles, SharedArrays, Random, LinearAlgebra, Clustering, Base
using DataFrames, CSV
using Statistics, ArgParse

nmats = 100 # iteraciones para samplear los eigenvalues de las correlaciones
n_kmedoids = 100 # iteraciones para calcular el minimo de kmedoids
n_kmedoids_out = 10 # Output de de las iteraciones del minimo de los kmedoids

function max_eigv(gene) #returns largest eigenvalue of a randomized matrix generated from data.
    gran = zeros(size(gene))
    for i = 1:size(gene,2)
        gran[:,i] = gene[shuffle(1:end),i] #Shuffle over genes
    end
 #    @time mr = corr_mat(gran)
    # println(mr[[end-1],end])
    mr = cor(transpose(gran))
    # println(mr[[end-1],end])
    return maximum(eigvals(mr))
end

#############
## MAIN
#############

## JULIA_NUM_THREADS=$(nproc)
# julia Julia/02-clusters-chrom.jl -c Data/clean/Healthy-clean.tsv

s = ArgParseSettings()
@add_arg_table! s begin
    # "-c"
    #     help = "Flag to calculate number of clusters"
    #     action = :store_true
    "expfile"
        help = "Expression Gene Matrix (read a clean format exp matrix)"
        arg_type = String
        required = true
end
parsed_args = parse_args(s)
# parsed_args = Dict([("c", true), ("expfile", "Data/clean/Basal-clean.tsv")])
println(parsed_args)
fname = basename(parsed_args["expfile"])
subtype = split(fname,'-')[1]
println("Working on " * subtype)
flush(stdout)

Random.seed!(34568)

##============ 0.5 - Input
#! gname s1 s2 ... gstart chromname
df = CSV.read(parsed_args["expfile"], delim='\t')

#! Selection per chromosome
for g in groupby(df, :chromname)
    cname = unique(g.chromname)[1]
    println("Working on chromosome: ", cname)
    #! Extract expression by chromosome
    gene = convert(Array{Float64,2},g[:, 2:end-2])
    println(size(gene))

    ##============ 01 - First correlation
    println("First correlation (m)")
    @time  m = cor(transpose(gene))

    ##============ 02 - Eigenvalues
    println("Get Eigenvalues from ",nmats," shuffles")
    l_evmax = SharedArray{Float64, 1}((nmats))
    @time for i in 1:nmats #doing the surrogate matrices and saving max eigv.
        # @time l_evmax[i] = max_eigv(gene)
        l_evmax[i] = max_eigv(gene)
        # println("Iteration: ",i,", max eigenvalue: ",l_evmax[i])
        flush(stdout)
    end
    #! Get eigenvalue of gene
    m_ev = filter(x->x>0.0000001, eigvals(m)) # Drop Zeros
    #! Todos los eigenvalues de la correlacion que sean mayores al mayor de las
    #! correlaciones muestreadas; el numero de estos eigenvalues es el numero de cluesters
    #! nc = 13, for Basal-All.txt with nmats=100
    nc = length(filter(x->x>maximum(l_evmax),m_ev))
    println("Number of clusters (nc): ",nc)

    ##============ 03 - KMedoids
    println("Starting kmedoids sampling")
    #! Make m in [0,1] insted [-1,1]
    dmat = map(x-> 1-abs(x),m)
    print(size(dmat))
    ens_c = []
    for i = 1:n_kmedoids
        push!(ens_c, kmedoids(dmat,nc))
        if (i % n_kmedoids_out) == 0
          println(i)
        end
    end

    #! De todas los clustereos cual es el que tiene el costo total minimo
    inx_c = findmin(map(x -> x.totalcost, ens_c))[2]
    c = ens_c[inx_c]
    println("Cluster counts: ",c.counts)

    ##============ 04 - Saving
    df = DataFrame([g[:,1],c.assignments,c.costs,g[:,end-1],g[:,end]])
    header = ["gname","clusterid","assignmentcost","gstart","chromosome"]
    rename!(df, Symbol.(header))
    odir = "Data/Clustered/" * subtype * "/"
    mkpath(odir)
    fname = odir*subtype*"-chr"*cname*"-clusters.tsv"
    CSV.write(fname,df)
end

max_stars_repo_name: josemaz/gene-matrices

max_issues_repo_name: josemaz/gene-matrices

max_forks_repo_name: josemaz/gene-matrices


text:
using Test, DWDataReader

const dir = joinpath(dirname(pathof(DWDataReader)), "..", "test", "testfiles")

@testset "DWDataReader" begin

    DWInit()
    @test DWDataReader.DWGetVersion() == 4020020
    DWDeInit()

    @testset "DWDataReader.File" begin

        include("basics.jl")
        include("testfiles.jl")

    end # @testset "DWDataReader.File"

end


text:
include("src/julia/string_bits_loading/loader_fns.jl")
xvar = [randstring(rand(1:8)) for i in 1:100_000];
x = rand(xvar, 10_000_000);
using BenchmarkTools
# @btime leading_bitsdot($x);
# @btime leading_bits_with_fast_pathdot($x);
@btime unsafe_loadbitsdot($x);          # 70ms
@btime load_bits_with_paddingdot($x);   # 240ms
# @btime copy_bits_fast_pathdot($x)     # can't test
@btime all_lensdot($x);                 # 190ms
@btime check_boundarydot($x);           # 72ms

xfixed = [randstring(8) for i in 1:100_000];
x = rand(xfixed, 10_000_000);
using BenchmarkTools
# @btime leading_bitsdot($x);
# @btime leading_bits_with_fast_pathdot($x);
@btime unsafe_loadbitsdot($x);          # 67ms
@btime load_bits_with_paddingdot($x);   # 78ms
# @btime copy_bits_fast_pathdot($x)     # can't test
@btime all_lensdot($x);                 # 81ms
@btime check_boundarydot($x);           # 74ms

##########################################################################
# Larger tests
######################

xvar = [randstring(rand(1:8)) for i in 1:1_000_000];
x = rand(xvar, 100_000_000);
using BenchmarkTools
# @btime leading_bitsdot($x);
# @btime leading_bits_with_fast_pathdot($x);
@btime unsafe_loadbitsdot($x);          # 70ms
# @btime load_bits_with_paddingdot($x);   # 240ms
# @btime copy_bits_fast_pathdot($x)     # can't test
# @btime all_lensdot($x);                 # 190ms
@btime check_boundarydot($x);           # 72ms


text:
# Spherical and cartesian coordinate conversions

"""
    azimuth(lon1, lat1, lon2, lat2, degrees::Bool=true) -> az

Compute the azimuth `az` from point (`lon1`, `lat1`) to (`lon2`, `lat2`) on the sphere.
Points and azimuth are read and returned in degrees by default; use `degrees=false`
for radians.
"""
function azimuth(lon1, lat1, lon2, lat2, degrees::Bool=true)
    if degrees
        lon1, lat1, lon2, lat2 = deg2rad(lon1), deg2rad(lat1), deg2rad(lon2), deg2rad(lat2)
    end
    azimuth = atan(sin(lon2-lon1)*cos(lat2),
                   cos(lat1)*sin(lat2) - sin(lat1)*cos(lat2)*cos(lon2-lon1))
    degrees ? rad2deg(azimuth) : azimuth
end

"""
    cart2geog(x, y, z, degrees::Bool=true) -> lon, lat, r

Compute the longitude, latitude and radius given the cartesian coordinates `x`,
`y` and `z`, where `x` is at (lon,lat) = (0,0), `y` is at (90°,0) and `z` is
through lat = 90°.
"""
function cart2geog(x, y, z, degrees::Bool=true)
    r = sqrt(x^2 + y^2 + z^2)
    r == 0. && return zero(x), zero(x), zero(x)
    lon = atan(y, x)
    lat = asin(z/r)
    degrees ? (rad2deg(lon), rad2deg(lat), r) : (lon, lat, r)
end

function cart2geog(x::AbstractArray, y::AbstractArray, z::AbstractArray, degrees::Bool=true)
    dims = size(x)
    dims == size(y) == size(z) || throw(ArgumentError("All arrays must have same length"))
    T = promote_type(float.(eltype.((x, y, z)))...)
    lon, lat, r = Array{T}(dims), Array{T}(dims), Array{T}(dims)
    for i in eachindex(lon)
        lon[i], lat[i], r[i] = cart2geog(x[i], y[i], z[i], degrees)
    end
    lon, lat, r
end

"""
    delta(lon1, lat1, lon2, lat2, degrees::Bool=true) -> d

Compute the angular distance `d` on the sphere between two points, (`lon1`, `lat1`)
and (`lon2`, `lat2`).  Points and distance are read and returned in degrees
by default; use `degrees=false` for radians.
"""
function delta(lon1, lat1, lon2, lat2, degrees::Bool=true)
    if degrees
        lon1, lat1, lon2, lat2 = deg2rad(lon1), deg2rad(lat1), deg2rad(lon2), deg2rad(lat2)
    end
    d = atan(sqrt(
               (cos(lat2)*sin(lon2-lon1))^2 + (cos(lat1)*sin(lat2) -
                sin(lat1)*cos(lat2)*cos(lon2-lon1))^2),
               sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon2-lon1)
              )
    degrees ? rad2deg(d) : d
end

"""
    geog2cart(lon, lat, r, degrees::Bool=true) -> x, y, z
    geog2cart(lon, lat, degrees::Bool=true) -> x, y, z

Return the cartesian coordinates given the geographic longitude, latitude and
radius `lon`, `lat` and `r`.

If `r` is not given, points are returned on the unit sphere.
"""
function geog2cart(lon, lat, r, degrees::Bool=true)
    points_valid(lon, lat, degrees) || error("geog2cart: Points are not on the sphere")
    degrees && begin lon, lat = deg2rad(lon), deg2rad(lat) end
    x = r*cos(lon)*cos(lat)
    y = r*sin(lon)*cos(lat)
    z = r*sin(lat)
    x, y, z
end
function geog2cart(lon::AbstractArray, lat::AbstractArray, r::AbstractArray, degrees::Bool=true)
    size(lon) == size(lat) == size(r) ||
        throw(ArgumentError("Sizes of lon, lat and r must be the same"))
    T = promote_type(float.(eltype.((lon, lat, r)))...)
    x = Array{T}(undef, size(lon))
    y, z = similar(x), similar(x)
    for i in eachindex(lon)
        x[i], y[i], z[i] = geog2cart(lon[i], lat[i], r[i], degrees)
    end
    x, y, z
end
geog2cart(lon::AbstractArray, lat::AbstractArray, degrees::Bool=true) =
    geog2cart(lon, lat, fill(one(eltype(lon)), size(lon)), degrees)

"""
    points_valid(lon, lat, degrees::Bool=true) -> ::Bool

Return `true` if all points in arrays `lon` and `lat` are on the sphere.
`lon` is ignored, but `lat` is checked to see if points are in the range
-90°–90° (-π–π).  Points are read and returned in degrees by default; use
`degrees=false` for radians.
"""
points_valid(lon, lat, degrees::Bool=true) =
    degrees ? !any(abs.(lat) .> 90.) : !any(abs.(lat) .> pi/2.)

"""
    geodesic_endpoint(lon, lat, az, delta, degrees::Bool=true) -> lon1, lat1

Compute the end point (`lon1`, `lat1`) reached by travelling on the sphere along
azimuth `az` for `delta` angular distance.  Points, angles and distance are read
and returned in degrees by default; use `degrees=false` for radians.
"""
function geodesic_endpoint(lon, lat, az, delta, degrees::Bool=true)
    if degrees
        lon, lat, az, delta = deg2rad(lon), deg2rad(lat), deg2rad(az), deg2rad(delta)
    end
    lat2 = asin(sin(lat)*cos(delta) + cos(lat)*sin(delta)*cos(az))
    lon2 = lon + atan(sin(az)*sin(delta)*cos(lat),
                      cos(delta)-sin(lat)*sin(lat2))
    if degrees
        lon2, lat2 = rad2deg(lon2), rad2deg(lat2)
    end
    lon2, lat2
end

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
"""

"""

struct ISolver <: AbstractSolver
	options::Options

	function ISolver(op::Options)
		new(op)
	end
end

ISolver(op::Pair{Symbol,<:Any}...) =  ISolver(Options(op))

init(S::AbstractSystem, O::Options) = init!(S, O)

function init!(S::AbstractSystem, O::Options)
	solver = ISolver(O)
	return solver
end

function preprocess(S::AbstractSystem, O::Options)
	type = eltype(S.X)

	dim_x = state_dim(S)

	solver = init(S, O)

	max_step = solver.options.dict[:max_step]::Int
	#TODO implement inerval arithmetic to avoid sampling
	Xsp = sampleX(solver.options.dict[:Xsampletype],solver.options.dict[:nXsamples],dim_x, type)

	if has_input(S)
		dim_u = input_dim(S)
		Usp = sampleUW(solver.options.dict[:nUsamples], dim_u, low(S.U), high(S.U))
	else
		dim_u = 1
		Usp = [[0.0]]
	end

	if has_uncertainty(S)
		dim_w = uncertainty_dim(S)
		Wsp = sampleUW(solver.options.dict[:nWsamples], dim_w, low(S.W), high(S.W))
	else
		dim_w = 1
		Wsp = [[0.0]]
	end

	iset = split(state_set(S),ones(Int,dim_x))

	return dim_x, max_step, Xsp, Usp, Wsp, iset
end

function computeISet(S::AbstractSystem, O::Options)
	dim_x, max_step, Xsp, Usp, Wsp, iset = preprocess(S, O)
	i = 0
	@inbounds while i < max_step
		println("Iteration: $i")
		number_of_cells = length(iset)
		println("Number of cells: $number_of_cells")
		println("")

		# select all cells to divide
		to_divide = collect(1:length(iset))
		# divide all cells
		subdivide_cells!(iset, to_divide, dim_x, i%dim_x + 1)
		# uncertainty set limit
		@inbounds for w in Wsp
			# create symbolic image (digraph)
			G = create_symbolic_image(iset, S, Xsp, Usp, w)
			# select nonleaving cells
			nlc_w = analyze_graph(G)
			# find intersection
			nlc = intersect(collect(1:length(iset)),nlc_w)
			# select nonleaving cells and discard leaving cells
			iset = iset[nlc]
		end
		# increment i
		i += 1
	end
	return Solution(iset, O)
end


text:
using LatentDiffEq
using FileIO
using Parameters: @with_kw
using ProgressMeter
using Random
using Statistics
using MLDataUtils
using BSON: @save
using Flux.Data: DataLoader
using Flux
using DiffEqFlux
using OrdinaryDiffEq
using StochasticDiffEq
using DiffEqSensitivity
using ModelingToolkit
using Images
using Plots
import GR

include("pendulum.jl")
include("nODE.jl")
include("create_data.jl")

################################################################################
## Arguments for the train function
@with_kw mutable struct Args
    ## Global model
    # model_type = GOKU()
    model_type = LatentODE()

    ## Latent Differential Equations
    # diffeq = Pendulum_friction()
    # diffeq = Pendulum()
    # diffeq = SPendulum()
    # diffeq = Pendulum_NN_friction()
    # diffeq = NODE(2, augment_dim = 2)
    diffeq = NODE(16)

    ## Training params
    η = 1e-3                        # learning rate
    decay = 0.0001f0                # decay applied to weights during optimisation
    batch_size = 64                 # minibatch size
    seq_len = 50                    # sequence length for training samples
    epochs = 1500                   # number of epochs for training
    seed = 1                        # random seed
    cuda = false                    # GPU usage (not working well yet)
    dt = 0.05                       # timestep for ode solve
    variational = true              # variational or deterministic training

    ## Annealing schedule
    start_β = 0f0                   # start value
    end_β = 1f0                     # end value
    n_cycle = 4                     # number of annealing cycles
    ratio = 0.9                     # proportion used to increase β (and 1-ratio used to fix β)

    ## Progressive observation training
    progressive_training = false    # progressive training usage
    prog_training_duration = 200    # number of epochs to reach the final seq_len
    start_seq_len = 5               # training sequence length at first step

    ## Visualization
    vis_len = 60                    # number of frames to visualize after each epoch
    save_figure = false             # true: save visualization figure in save_path folder
                                    # false: display image instead of saving it
end

################################################################################
################################################################################
## Training done manualy

function train(; kws...)
    ## Load hyperparameters and GPU config
    args = Args(; kws...)
    @unpack_Args args

    seed > 0 && Random.seed!(seed)

    device = cpu
    @info "Training on CPU"

    ############################################################################
    ## Prepare training data

    root_dir = @__DIR__
    data_path = "$root_dir/data/data.bson"

    if ~isfile(data_path)
        @info "Generating data"
        latent_data, u0s, ps, high_dim_data = generate_dataset(diffeq = diffeq)
        data = (latent_data, u0s, ps, high_dim_data)
        mkpath("$root_dir/data")
        @save data_path data
    end

    data_loaded = load(data_path, :data)
    train_data = data_loaded[4]
    latent_data = data_loaded[1]

    # stack time for each sample
    train_data = Flux.stack.(train_data, 3)

    # stack all samples
    train_data = Flux.stack(train_data, 4) # 28x28x400x450
    h, w, full_seq_len, observations = size(train_data)
    latent_data = Flux.stack(latent_data, 3)

    # vectorize frames
    train_data = reshape(train_data, :, full_seq_len, observations) # input_dim, time_size, samples
    train_data = Float32.(train_data)

    train_set, val_set = Array.(splitobs(train_data, 0.9))
    train_set_latent, val_set_latent = Array.(splitobs(latent_data, 0.9))

    # loader_train = DataLoader(Array(train_set), batchsize=batch_size, shuffle=true, partial=false)
    loader_train = DataLoader((train_set, train_set_latent), batchsize=batch_size, shuffle=true, partial=false)
    val_set_time_unstacked = Flux.unstack(val_set, 2)

    input_dim = size(train_set,1)

    ############################################################################
    # Create model

    encoder_layers, decoder_layers = default_layers(model_type, input_dim, diffeq, device = device)
    model = LatentDiffEqModel(model_type, encoder_layers, decoder_layers)

    # Get parameters
    ps = Flux.params(model)

    ############################################################################
    ## Define optimizer
    # opt = AdaBelief(η)
    # opt = ADAM(η)
    opt = ADAMW(η,(0.9,0.999), decay)

    ############################################################################
    ## Various definitions
    if progressive_training
        prog_seq_lengths = range(start_seq_len, seq_len, step=(seq_len-start_seq_len)/(prog_training_duration-1))
        prog_seq_lengths = Int.(round.(prog_seq_lengths))
    else
        prog_training_duration = 0
    end

    # KL annealing scheduling
    annealing_schedule = frange_cycle_linear(epochs, start_β, end_β, n_cycle, ratio)

    # Preparation for saving best models weights
    mkpath("$root_dir/output")
    best_val_loss = Inf32
    val_loss = 0f0

    # mkpath("$root_dir/output")
    # args = struct2dict(args)
    # @save "$root_dir/output/args.bson" args

    ## Visualization options
    if save_figure
        mkpath("$root_dir/output/visualization")
        GR.inline("pdf")
    end
    ############################################################################
    ## Main train loop
    @info "Start Training of $(typeof(model_type))-net, total $epochs epochs"
    for epoch = 1:epochs

        # Set annealing factor
        β = annealing_schedule[epoch]

        ## set a sequence length for training samples
        seq_len = epoch ≤ prog_training_duration ? prog_seq_lengths[epoch] : seq_len

        # Model evaluation length
        t = range(0.f0, step=dt, length=seq_len)

        @info "Epoch $epoch .. (Sequence training length $seq_len)"
        progress = Progress(length(loader_train))

        for data in loader_train
            x, latent = data

            # Use only random sequences of length seq_len for the current minibatch
            x = time_loader(x, full_seq_len, seq_len)

            # Run the model with the current parameters and compute the loss
            loss, back = Flux.pullback(ps) do
                loss_batch(model, x |> device, t, β, variational)
            end

            # Backpropagate and update
            grad = back(1f0)
            Flux.Optimise.update!(opt, ps, grad)

            # Use validation set to get loss and visualisation
            t_val = range(0.f0, step=dt, length=length(val_set_time_unstacked))
            val_loss = loss_batch(model, val_set_time_unstacked |> device, t_val, β, false)

            # progress meter
            next!(progress; showvalues=[(:loss, loss),(:val_loss, val_loss)])
        end

        if device != gpu
            visualize_val_image(model, val_set |> device, val_set_latent, vis_len, dt, h, w, save_figure)
        end

        if val_loss < best_val_loss
            best_val_loss = deepcopy(val_loss)
            weights = Flux.params(model)
            @save "$root_dir/output/best_model_weights.bson" weights
            @info "Model saved"
        end
    end
end

################################################################################
## Loss definition
# f(x) = 0.0000000001f0/(x^3)
function loss_batch(model, x, t, β, variational)

    # Make prediction
    X̂, μ, logσ² = model(x, t, variational)
    x̂, ẑ, l̂ = X̂

    # Compute reconstruction loss
    reconstruction_loss = vector_mse(x, x̂)

    # # Compute 1st order difference loss
    # x_diff1 = diff(x)
    # x̂_diff1 = diff(x̂)
    # diff1_loss = vector_mse(x_diff1, x̂_diff1)
    # # Compute 2st order difference loss
    # x_diff2 = diff(x_diff1)
    # x̂_diff2 = diff(x̂_diff1)
    # diff2_loss = vector_mse(x_diff2, x̂_diff2)

    # Compute KL losses from parameter and initial value estimation
    kl_loss = vector_kl(μ, logσ²)

    return reconstruction_loss + β * kl_loss #+ λ*diff1_loss + λ*diff2_loss
end

################################################################################
## Visualization function

function visualize_val_image(model, val_set, val_set_latent, vis_len, dt, h, w, save_figure)
    j = rand(1:size(val_set,3))
    idxs = rand_time(size(val_set,2), vis_len)
    X_test = val_set[:, idxs, j]
    true_latent = val_set_latent[:,idxs,j]

    frames_test = [Gray.(reshape(x,h,w)) for x in eachcol(X_test)]
    X_test = reshape(X_test, Val(3))
    x = Flux.unstack(X_test, 2)
    t_val = range(0.f0, step=dt, length=vis_len)

    X̂, μ, logσ² = model(x, t_val)
    x̂, ẑ, l̂ = X̂

    ẑ = Flux.stack(ẑ, 2)

    plt1 = plot(ẑ[:,:,1]', label="inferred",legend=:topleft, rightmargin = 2.0Plots.cm)
    ylabel!("State variables")
    xlabel!("time")
    plt1 = plot!(twinx(), true_latent[1,:], color=:green, box = :on, xticks=:none, label="ground truth")

    x̂ = Flux.stack(x̂, 2)
    frames_pred = [Gray.(reshape(x,h,w)) for x in eachslice(x̂, dims=2)]

    frames_test = frames_test[1:6:end]
    frames_pred = frames_pred[1:6:end]

    plt2 = mosaicview(frames_test..., frames_pred..., nrow=2, rowmajor=true)
    # plt2 = plot(plt2, leg = false, ticks = nothing, border = :none)
    # plt = plot(plt1, plt2, layout = @layout([a; b]))
    # save_figure ? savefig(plt, "output/visualization/fig.pdf") : display(plt)

    display(plt1)
    display(plt2)
end

train()


text:
using BenchmarkTools
using FoldsCatwalk
using Transducers

returnnothing(_) = nothing
returnnothing(_, _) = nothing

foreach1(f, xs, exc) = transduce(Map(f), returnnothing, nothing, xs, exc)

type_instability(x) = Val(round(Int, 2x))

trydispatch(_) = nothing
trydispatch(::Val{typemax(Int)}) = error("unexpected!") # so that Julia has to try dispatch

function foreach_baseline(xs)
    itr = xs |> Map(type_instability)
    foreach1(trydispatch, itr, SequentialEx())
end

function foreach_catwalk(xs)
    itr = xs |> Map(type_instability) |> OptimizeInner()
    foreach1(trydispatch, itr, CatwalkEx())
end

function demo_foreach()
    xs = randn(1000_000)
    @btime foreach_baseline($xs)
    @btime foreach_catwalk($xs)
end


text:
"""
        fs_expansion(x, u, v, L, T, n = i -> 2^-(1 + i/2))
Global recursion: find the value of the process on the dyadic points given by the first L Levels,
`x` is a vector containing coordinates `u`,`v` initial and final point
`T` length of the basis
`n` is the normalization function (in our case known)
"""
function fs_expansion(x::Vector{Float64}, u::Float64, v::Float64, L::Int64, T::Float64, n = i -> 2^-(1 + i/2))
        y = fill(NaN, 2<<L + 1)
        y[1] = u
        y[end] = v
        k = 1
        c = sqrt(T)
        for i in 0:L
                for j in 0:(1<<i - 1)
                        l = 1<<(L-i) + 2<<(L-i)*(j) + 1 #new index
                        y[l] = x[k]*n(i)*c + (y[l + 1<<(L-i)] + y[l - 1<<(L-i)])/2
                        k += 1
                end
        end
        y
end

"""
        fs_expansion(ϕ::Fs, ξ, u, v, L::Int64, T::Float64, n = i -> 2^-(1 + i/2))
Local recursion, like Global but on the support of the Faber Schauder function `ϕ` up to level `L`
This function evaluates only the necessary functions for the evaluation of the points inside the support
of the function ϕ.
`x` is a vector containing coordinates `u`,`v` initial and final point
`T` length of the basis
`n` is the normalization function (in our case known)

"""
function fs_expansion(ϕ::Fs, ξ::Vector{Float64}, u::Float64, v::Float64, L::Int64, T::Float64, n = i -> 2^-(1 + i/2))
        y = fill(NaN, 2<<(L) + 1)   #number of points in which we require the evaluation of the function
        y[end] = v
        y[1] = u
        c = sqrt(T)
        for k in ϕ.nhb
                i,j = Faber(k)
                l = 2<<(L-i-1) + 2<<(L-i)*(j) + 1 #new index
                y[l] = ξ[k]*n(i)*c + (y[l + 2<<(L-i-1)]  + y[l - 2<<(L-i-1)])/2
                k += 1
        end
        y[(2<<(L-ϕ.i)*ϕ.j + 1):(2<<(L-ϕ.i)*(ϕ.j + 1) + 1)]
end

"""
        interpolate(x:: array, y::array, x0)

simple interpolation function given `x`, `y`: 2-element array, `x0` floating
"""
function interpolate(x::Vector{Float64}, y::Vector{Float64}, x0::Float64)
        if x0 == x[1]
                return y[1]
        elseif x0 == x[2]
                return y[2]
        else
                return ( y[2] - (y[2] - y[1])/(x[2] - x[1])*(x[2] - x0))
        end
end

"""
        fs_expansion(S::System, t::Float64, u, v, n = i -> 2^-(1 + i/2))

Find value of the process (Piecewise linear) for any t ∈ [0, T]
does not evaluate the  whole path, but just the points needed for t
"""
function fs_expansion(S::System, t::Float64, u::Float64, v::Float64, n = i -> 2^-(1 + i/2))
        dt = 0:S.T/(2<<S.L):S.T
        k = (searchsortedfirst(dt, t) - 1)
        j0 =  Int(ceil(k/2))-1
        n0 = Faber(S.L, j0)
        if k % 2 != 0
                return interpolate([dt[k], dt[k + 1]], fs_expansion(S.ϕ[n0], S.ξ, u, v, S.L, S.T)[1:2], t)
        else
                return interpolate([dt[k], dt[k + 1]], fs_expansion(S.ϕ[n0], S.ξ, u, v,  S.L, S.T)[2:3], t)
        end
end

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
const v0 = -84.624
const v1 = 10.0
const C_K1 = 1.0f0
const C_x1 = 1.0f0
const C_Na = 1.0f0
const C_s = 1.0f0
const D_Ca = 0.0f0
const D_Na = 0.0f0
const g_s = 0.09f0
const g_Na = 4.0f0
const g_NaC = 0.005f0
const ENa = 50.0f0 + D_Na
const γ = 0.5f0
const C_m = 1.0f0

mutable struct BeelerReuterCpu <: Function
    t::Float64              # the last timestep time to calculate Δt
    diff_coef::Float64      # the diffusion-coefficient (coupling strength)

    C::Array{Float32, 2}    # intracellular calcium concentration
    M::Array{Float32, 2}    # sodium current activation gate (m)
    H::Array{Float32, 2}    # sodium current inactivation gate (h)
    J::Array{Float32, 2}    # sodium current slow inactivaiton gate (j)
    D::Array{Float32, 2}    # calcium current activaiton gate (d)
    F::Array{Float32, 2}    # calcium current inactivation gate (f)
    XI::Array{Float32, 2}   # inward-rectifying potassium current (iK1)

    Δu::Array{Float64, 2}   # place-holder for the Laplacian

    function BeelerReuterCpu(u0, diff_coef)
        self = new()

        ny, nx = size(u0)
        self.t = 0.0
        self.diff_coef = diff_coef

        self.C = fill(0.0001f0, (ny,nx))
        self.M = fill(0.01f0, (ny,nx))
        self.H = fill(0.988f0, (ny,nx))
        self.J = fill(0.975f0, (ny,nx))
        self.D = fill(0.003f0, (ny,nx))
        self.F = fill(0.994f0, (ny,nx))
        self.XI = fill(0.0001f0, (ny,nx))

        self.Δu = zeros(ny,nx)

        return self
    end
end

# 5-point stencil
function laplacian(Δu, u)
    n1, n2 = size(u)

    # internal nodes
    for j = 2:n2-1
        for i = 2:n1-1
            @inbounds  Δu[i,j] = u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i,j-1] - 4*u[i,j]
        end
    end

    # left/right edges
    for i = 2:n1-1
        @inbounds Δu[i,1] = u[i+1,1] + u[i-1,1] + 2*u[i,2] - 4*u[i,1]
        @inbounds Δu[i,n2] = u[i+1,n2] + u[i-1,n2] + 2*u[i,n2-1] - 4*u[i,n2]
    end

    # top/bottom edges
    for j = 2:n2-1
        @inbounds Δu[1,j] = u[1,j+1] + u[1,j-1] + 2*u[2,j] - 4*u[1,j]
        @inbounds Δu[n1,j] = u[n1,j+1] + u[n1,j-1] + 2*u[n1-1,j] - 4*u[n1,j]
    end

    # corners
    @inbounds Δu[1,1] = 2*(u[2,1] + u[1,2]) - 4*u[1,1]
    @inbounds Δu[n1,1] = 2*(u[n1-1,1] + u[n1,2]) - 4*u[n1,1]
    @inbounds Δu[1,n2] = 2*(u[2,n2] + u[1,n2-1]) - 4*u[1,n2]
    @inbounds Δu[n1,n2] = 2*(u[n1-1,n2] + u[n1,n2-1]) - 4*u[n1,n2]
end

@inline function rush_larsen(g, α, β, Δt)
    inf = α/(α+β)
    τ = 1f0 / (α+β)
    return clamp(g + (g - inf) * expm1(-Δt/τ), 0f0, 1f0)
end

function update_M_cpu(g, v, Δt)
    # the condition is needed here to prevent NaN when v == 47.0
    α = isapprox(v, 47.0f0) ? 10.0f0 : -(v+47.0f0) / (exp(-0.1f0*(v+47.0f0)) - 1.0f0)
    β = (40.0f0 * exp(-0.056f0*(v+72.0f0)))
    return rush_larsen(g, α, β, Δt)
end

function update_H_cpu(g, v, Δt)
    α = 0.126f0 * exp(-0.25f0*(v+77.0f0))
    β = 1.7f0 / (exp(-0.082f0*(v+22.5f0)) + 1.0f0)
   return rush_larsen(g, α, β, Δt)
end

function update_J_cpu(g, v, Δt)
    α = (0.55f0 * exp(-0.25f0*(v+78.0f0))) / (exp(-0.2f0*(v+78.0f0)) + 1.0f0)
    β = 0.3f0 / (exp(-0.1f0*(v+32.0f0)) + 1.0f0)
    return rush_larsen(g, α, β, Δt)
end

function update_D_cpu(g, v, Δt)
    α = γ * (0.095f0 * exp(-0.01f0*(v-5.0f0))) / (exp(-0.072f0*(v-5.0f0)) + 1.0f0)
    β = γ * (0.07f0 * exp(-0.017f0*(v+44.0f0))) / (exp(0.05f0*(v+44.0f0)) + 1.0f0)
    return rush_larsen(g, α, β, Δt)
end

function update_F_cpu(g, v, Δt)
    α = γ * (0.012f0 * exp(-0.008f0*(v+28.0f0))) / (exp(0.15f0*(v+28.0f0)) + 1.0f0)
    β = γ * (0.0065f0 * exp(-0.02f0*(v+30.0f0))) / (exp(-0.2f0*(v+30.0f0)) + 1.0f0)
    return rush_larsen(g, α, β, Δt)
end

function update_XI_cpu(g, v, Δt)
    α = (0.0005f0 * exp(0.083f0*(v+50.0f0))) / (exp(0.057f0*(v+50.0f0)) + 1.0f0)
    β = (0.0013f0 * exp(-0.06f0*(v+20.0f0))) / (exp(-0.04f0*(v+20.0f0)) + 1.0f0)
    return rush_larsen(g, α, β, Δt)
end

function update_C_cpu(g, d, f, v, Δt)
    ECa = D_Ca - 82.3f0 - 13.0278f0 * log(g)
    kCa = C_s * g_s * d * f
    iCa = kCa * (v - ECa)
    inf = 1.0f-7 * (0.07f0 - g)
    τ = 1f0 / 0.07f0
    return g + (g - inf) * expm1(-Δt/τ)
end

function update_gates_cpu(u, XI, M, H, J, D, F, C, Δt)
    let Δt = Float32(Δt)
        n1, n2 = size(u)
        for j = 1:n2
            for i = 1:n1
                v = Float32(u[i,j])

                XI[i,j] = update_XI_cpu(XI[i,j], v, Δt)
                M[i,j] = update_M_cpu(M[i,j], v, Δt)
                H[i,j] = update_H_cpu(H[i,j], v, Δt)
                J[i,j] = update_J_cpu(J[i,j], v, Δt)
                D[i,j] = update_D_cpu(D[i,j], v, Δt)
                F[i,j] = update_F_cpu(F[i,j], v, Δt)

                C[i,j] = update_C_cpu(C[i,j], D[i,j], F[i,j], v, Δt)
            end
        end
    end
end

# iK1 is the inward-rectifying potassium current
function calc_iK1(v)
    ea = exp(0.04f0*(v+85f0))
    eb = exp(0.08f0*(v+53f0))
    ec = exp(0.04f0*(v+53f0))
    ed = exp(-0.04f0*(v+23f0))
    return 0.35f0 * (4f0*(ea-1f0)/(eb + ec)
            + 0.2f0 * (isapprox(v, -23f0) ? 25f0 : (v+23f0) / (1f0-ed)))
end

# ix1 is the time-independent background potassium current
function calc_ix1(v, xi)
    ea = exp(0.04f0*(v+77f0))
    eb = exp(0.04f0*(v+35f0))
    return xi * 0.8f0 * (ea-1f0) / eb
end

# iNa is the sodium current (similar to the classic Hodgkin-Huxley model)
function calc_iNa(v, m, h, j)
    return C_Na * (g_Na * m^3 * h * j + g_NaC) * (v - ENa)
end

# iCa is the calcium current
function calc_iCa(v, d, f, c)
    ECa = D_Ca - 82.3f0 - 13.0278f0 * log(c)    # ECa is the calcium reversal potential
    return C_s * g_s * d * f * (v - ECa)
end

function update_du_cpu(du, u, XI, M, H, J, D, F, C)
    n1, n2 = size(u)

    for j = 1:n2
        for i = 1:n1
            v = Float32(u[i,j])

            # calculating individual currents
            iK1 = calc_iK1(v)
            ix1 = calc_ix1(v, XI[i,j])
            iNa = calc_iNa(v, M[i,j], H[i,j], J[i,j])
            iCa = calc_iCa(v, D[i,j], F[i,j], C[i,j])

            # total current
            I_sum = iK1 + ix1 + iNa + iCa

            # the reaction part of the reaction-diffusion equation
            du[i,j] = -I_sum / C_m
        end
    end
end

function (f::BeelerReuterCpu)(du, u, p, t)
    Δt = t - f.t

    if Δt != 0 || t == 0
        update_gates_cpu(u, f.XI, f.M, f.H, f.J, f.D, f.F, f.C, Δt)
        f.t = t
    end

    laplacian(f.Δu, u)

    # calculate the reaction portion
    update_du_cpu(du, u, f.XI, f.M, f.H, f.J, f.D, f.F, f.C)

    # ...add the diffusion portion
    du .+= f.diff_coef .* f.Δu
end

const N = 192;
u0 = fill(v0, (N, N));
u0[90:102,90:102] .= v1;   # a small square in the middle of the domain

using Plots
heatmap(u0)

using DifferentialEquations, Sundials

deriv_cpu = BeelerReuterCpu(u0, 1.0);
prob = ODEProblem(deriv_cpu, u0, (0.0, 50.0));

@time sol = solve(prob, CVODE_BDF(linear_solver=:GMRES), saveat=100.0);

heatmap(sol.u[end])

using CUDA

mutable struct BeelerReuterGpu <: Function
    t::Float64                  # the last timestep time to calculate Δt
    diff_coef::Float64          # the diffusion-coefficient (coupling strength)

    d_C::CuArray{Float32, 2}    # intracellular calcium concentration
    d_M::CuArray{Float32, 2}    # sodium current activation gate (m)
    d_H::CuArray{Float32, 2}    # sodium current inactivation gate (h)
    d_J::CuArray{Float32, 2}    # sodium current slow inactivaiton gate (j)
    d_D::CuArray{Float32, 2}    # calcium current activaiton gate (d)
    d_F::CuArray{Float32, 2}    # calcium current inactivation gate (f)
    d_XI::CuArray{Float32, 2}   # inward-rectifying potassium current (iK1)

    d_u::CuArray{Float64, 2}    # place-holder for u in the device memory
    d_du::CuArray{Float64, 2}   # place-holder for d_u in the device memory

    Δv::Array{Float64, 2}       # place-holder for voltage gradient

    function BeelerReuterGpu(u0, diff_coef)
        self = new()

        ny, nx = size(u0)
        @assert (nx % 16 == 0) && (ny % 16 == 0)
        self.t = 0.0
        self.diff_coef = diff_coef

        self.d_C = CuArray(fill(0.0001f0, (ny,nx)))
        self.d_M = CuArray(fill(0.01f0, (ny,nx)))
        self.d_H = CuArray(fill(0.988f0, (ny,nx)))
        self.d_J = CuArray(fill(0.975f0, (ny,nx)))
        self.d_D = CuArray(fill(0.003f0, (ny,nx)))
        self.d_F = CuArray(fill(0.994f0, (ny,nx)))
        self.d_XI = CuArray(fill(0.0001f0, (ny,nx)))

        self.d_u = CuArray(u0)
        self.d_du = CuArray(zeros(ny,nx))

        self.Δv = zeros(ny,nx)

        return self
    end
end

function rush_larsen_gpu(g, α, β, Δt)
    inf = α/(α+β)
    τ = 1.0/(α+β)
    return clamp(g + (g - inf) * CUDAnative.expm1(-Δt/τ), 0f0, 1f0)
end

function update_M_gpu(g, v, Δt)
    # the condition is needed here to prevent NaN when v == 47.0
    α = isapprox(v, 47.0f0) ? 10.0f0 : -(v+47.0f0) / (CUDAnative.exp(-0.1f0*(v+47.0f0)) - 1.0f0)
    β = (40.0f0 * CUDAnative.exp(-0.056f0*(v+72.0f0)))
    return rush_larsen_gpu(g, α, β, Δt)
end

function update_H_gpu(g, v, Δt)
    α = 0.126f0 * CUDAnative.exp(-0.25f0*(v+77.0f0))
    β = 1.7f0 / (CUDAnative.exp(-0.082f0*(v+22.5f0)) + 1.0f0)
    return rush_larsen_gpu(g, α, β, Δt)
end

function update_J_gpu(g, v, Δt)
    α = (0.55f0 * CUDAnative.exp(-0.25f0*(v+78.0f0))) / (CUDAnative.exp(-0.2f0*(v+78.0f0)) + 1.0f0)
    β = 0.3f0 / (CUDAnative.exp(-0.1f0*(v+32.0f0)) + 1.0f0)
    return rush_larsen_gpu(g, α, β, Δt)
end

function update_D_gpu(g, v, Δt)
    α = γ * (0.095f0 * CUDAnative.exp(-0.01f0*(v-5.0f0))) / (CUDAnative.exp(-0.072f0*(v-5.0f0)) + 1.0f0)
    β = γ * (0.07f0 * CUDAnative.exp(-0.017f0*(v+44.0f0))) / (CUDAnative.exp(0.05f0*(v+44.0f0)) + 1.0f0)
    return rush_larsen_gpu(g, α, β, Δt)
end

function update_F_gpu(g, v, Δt)
    α = γ * (0.012f0 * CUDAnative.exp(-0.008f0*(v+28.0f0))) / (CUDAnative.exp(0.15f0*(v+28.0f0)) + 1.0f0)
    β = γ * (0.0065f0 * CUDAnative.exp(-0.02f0*(v+30.0f0))) / (CUDAnative.exp(-0.2f0*(v+30.0f0)) + 1.0f0)
    return rush_larsen_gpu(g, α, β, Δt)
end

function update_XI_gpu(g, v, Δt)
    α = (0.0005f0 * CUDAnative.exp(0.083f0*(v+50.0f0))) / (CUDAnative.exp(0.057f0*(v+50.0f0)) + 1.0f0)
    β = (0.0013f0 * CUDAnative.exp(-0.06f0*(v+20.0f0))) / (CUDAnative.exp(-0.04f0*(v+20.0f0)) + 1.0f0)
    return rush_larsen_gpu(g, α, β, Δt)
end

function update_C_gpu(c, d, f, v, Δt)
    ECa = D_Ca - 82.3f0 - 13.0278f0 * CUDAnative.log(c)
    kCa = C_s * g_s * d * f
    iCa = kCa * (v - ECa)
    inf = 1.0f-7 * (0.07f0 - c)
    τ = 1f0 / 0.07f0
    return c + (c - inf) * CUDAnative.expm1(-Δt/τ)
end

# iK1 is the inward-rectifying potassium current
function calc_iK1(v)
    ea = CUDAnative.exp(0.04f0*(v+85f0))
    eb = CUDAnative.exp(0.08f0*(v+53f0))
    ec = CUDAnative.exp(0.04f0*(v+53f0))
    ed = CUDAnative.exp(-0.04f0*(v+23f0))
    return 0.35f0 * (4f0*(ea-1f0)/(eb + ec)
            + 0.2f0 * (isapprox(v, -23f0) ? 25f0 : (v+23f0) / (1f0-ed)))
end

# ix1 is the time-independent background potassium current
function calc_ix1(v, xi)
    ea = CUDAnative.exp(0.04f0*(v+77f0))
    eb = CUDAnative.exp(0.04f0*(v+35f0))
    return xi * 0.8f0 * (ea-1f0) / eb
end

# iNa is the sodium current (similar to the classic Hodgkin-Huxley model)
function calc_iNa(v, m, h, j)
    return C_Na * (g_Na * m^3 * h * j + g_NaC) * (v - ENa)
end

# iCa is the calcium current
function calc_iCa(v, d, f, c)
    ECa = D_Ca - 82.3f0 - 13.0278f0 * CUDAnative.log(c)    # ECa is the calcium reversal potential
    return C_s * g_s * d * f * (v - ECa)
end

function update_gates_gpu(u, XI, M, H, J, D, F, C, Δt)
    i = (blockIdx().x-UInt32(1)) * blockDim().x + threadIdx().x
    j = (blockIdx().y-UInt32(1)) * blockDim().y + threadIdx().y

    v = Float32(u[i,j])

    let Δt = Float32(Δt)
        XI[i,j] = update_XI_gpu(XI[i,j], v, Δt)
        M[i,j] = update_M_gpu(M[i,j], v, Δt)
        H[i,j] = update_H_gpu(H[i,j], v, Δt)
        J[i,j] = update_J_gpu(J[i,j], v, Δt)
        D[i,j] = update_D_gpu(D[i,j], v, Δt)
        F[i,j] = update_F_gpu(F[i,j], v, Δt)

        C[i,j] = update_C_gpu(C[i,j], D[i,j], F[i,j], v, Δt)
    end
    nothing
end

function update_du_gpu(du, u, XI, M, H, J, D, F, C)
    i = (blockIdx().x-UInt32(1)) * blockDim().x + threadIdx().x
    j = (blockIdx().y-UInt32(1)) * blockDim().y + threadIdx().y

    v = Float32(u[i,j])

    # calculating individual currents
    iK1 = calc_iK1(v)
    ix1 = calc_ix1(v, XI[i,j])
    iNa = calc_iNa(v, M[i,j], H[i,j], J[i,j])
    iCa = calc_iCa(v, D[i,j], F[i,j], C[i,j])

    # total current
    I_sum = iK1 + ix1 + iNa + iCa

    # the reaction part of the reaction-diffusion equation
    du[i,j] = -I_sum / C_m
    nothing
end

function (f::BeelerReuterGpu)(du, u, p, t)
    L = 16   # block size
    Δt = t - f.t
    copyto!(f.d_u, u)
    ny, nx = size(u)

    if Δt != 0 || t == 0
        @cuda blocks=(ny÷L,nx÷L) threads=(L,L) update_gates_gpu(
            f.d_u, f.d_XI, f.d_M, f.d_H, f.d_J, f.d_D, f.d_F, f.d_C, Δt)
        f.t = t
    end

    laplacian(f.Δv, u)

    # calculate the reaction portion
    @cuda blocks=(ny÷L,nx÷L) threads=(L,L) update_du_gpu(
        f.d_du, f.d_u, f.d_XI, f.d_M, f.d_H, f.d_J, f.d_D, f.d_F, f.d_C)

    copyto!(du, f.d_du)

    # ...add the diffusion portion
    du .+= f.diff_coef .* f.Δv
end

using DifferentialEquations, Sundials

deriv_gpu = BeelerReuterGpu(u0, 1.0);
prob = ODEProblem(deriv_gpu, u0, (0.0, 50.0));
@time sol = solve(prob, CVODE_BDF(linear_solver=:GMRES), saveat=100.0);

heatmap(sol.u[end])

using SciMLTutorials
SciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLTutorialsOutput

max_issues_repo_name: SciML/SciMLTutorialsOutput

max_forks_repo_name: SciML/SciMLTutorialsOutput


text:
function learn_optim(index::Int, model::Dict{String,Any}, training_df::DataFrame;
    pₒ::Union{Nothing,Array{Float64,1}} = nothing)

    # main simulation loop -
    SF = 1e9

    # setup static -
    sfa = model["static_factors_array"]
    sfa[1] = training_df[index,:TFPI]       # 1 TFPI
    sfa[2] = training_df[index,:AT]         # 2 AT
    sfa[3] = (5e-12) * SF               # 3 TF
    sfa[6] = 0.005                      # 6 TRAUMA

    # grab the multiplier from the data -
    ℳ = model["number_of_dynamic_states"]
    xₒ = zeros(ℳ)
    xₒ[1] = training_df[index, :II]         # 1 FII
    xₒ[2] = training_df[index, :VII]        # 2 FVII
    xₒ[3] = training_df[index, :V]          # 3 FV
    xₒ[4] = training_df[index, :X]          # 4 FX
    xₒ[5] = training_df[index, :VIII]       # 5 FVIII
    xₒ[6] = training_df[index, :IX]         # 6 FIX
    xₒ[7] = training_df[index, :XI]         # 7 FXI
    xₒ[8] = training_df[index, :XII]        # 8 FXII
    xₒ[9] = (1e-14)*SF                      # 9 FIIa
    model["initial_condition_vector"] = xₒ

    # what is the output array?
    Y = Array{Float64,1}(undef,5)
    Y[1] = training_df[index, :Lagtime]
    Y[2] = training_df[index, :Peak]
    Y[3] = training_df[index, Symbol("T.Peak")]
    Y[4] = training_df[index, :Max]
    Y[5] = training_df[index, :AUC]

    # setup initial parameter values and bounds array -
    κ = [

            # default: hand fit set -
            0.061   0.01 10.0   ; # 1
            1.0     0.01 10.0   ; # 2
            1.0     0.01 10.0   ; # 3
            1.0     0.01 10.0   ; # 4
            1.0     0.01 10.0   ; # 5
            1.0     0.01 10.0   ; # 6
            1.0     0.01 10.0   ; # 7
            1.0     0.01 10.0   ; # 8
            0.70    0.01 10.0   ; # 9
            0.11    0.01 10.0   ; # 10
            0.045   0.01 10.0   ; # 11
            0.065   0.01 10.0   ; # 12
        ];

    # set default set as the start -
    if (isnothing(pₒ) == true)
        P = length(κ[:,1])
        σ = 0.1 # move up to 10%
        pₒ = κ[:,1].*(1 .+ σ*rand(-1:1,P))
    end

    # setup the objective function -
    inner_optimizer = NelderMead()
    obj_function(p) =  loss_scalar(p, Y, model)
    results = optimize(obj_function, κ[:,2], κ[:,3], pₒ, Fminbox(inner_optimizer),
        Optim.Options(time_limit = 600, show_trace = true, show_every = 10, iterations=100))

    # grab the best parameters -
    p_best = Optim.minimizer(results)

    # run the sim w/the best parameters -
    # 1 - 9 : α vector
    model["α"] = p_best[1:9]
    G = model["G"]
    idx = indexin(model, "FVIIa")   # 10
    G[idx, 4] = p_best[10]
    idx = indexin(model, "AT")      # 11
    G[idx, 9] = p_best[11]
    idx = indexin(model, "TFPI")    # 12
    G[idx, 1] = -1*p_best[12]

    # run the model -
    (T,U) = evaluate(model)
    Xₘ = hcat(U...)
    Yₘ = model_output_vector(T, Xₘ[9,:]) # properties of the Thrombin curve

    return (p_best, T, Xₘ, Yₘ, Y)
end

function learn(index::Int; pₒ::Union{Nothing,Array{Float64,1}} = nothing)

    # let'e experiment with a simple evolutionary algorithm -
    NTRIALS = 12500
    σ = 0.05

    # load the training data -
    path_to_training_data = joinpath(_PATH_TO_DATA, "Training-Synthetic-Thrombin-TF-1K.csv")
    training_df = CSV.read(path_to_training_data, DataFrame)

    # build the model structure -
    path_to_model_file = joinpath(_PATH_TO_MODEL, "Coagulation.net")
    model_buffer = read_model_file(path_to_model_file)

    # build the default model structure -
    model = build_default_model_dictionary(model_buffer)

    # main simulation loop -
    SF = 1e9

    # setup static -
    sfa = model["static_factors_array"]
    sfa[1] = training_df[index,:TFPI]       # 1 TFPI
    sfa[2] = training_df[index,:AT]         # 2 AT
    sfa[3] = (5e-12) * SF               # 3 TF
    sfa[6] = 0.005                      # 6 TRAUMA

    # grab the multiplier from the data -
    ℳ = model["number_of_dynamic_states"]
    xₒ = zeros(ℳ)
    xₒ[1] = training_df[index, :II]         # 1 FII
    xₒ[2] = training_df[index, :VII]        # 2 FVII
    xₒ[3] = training_df[index, :V]          # 3 FV
    xₒ[4] = training_df[index, :X]          # 4 FX
    xₒ[5] = training_df[index, :VIII]       # 5 FVIII
    xₒ[6] = training_df[index, :IX]         # 6 FIX
    xₒ[7] = training_df[index, :XI]         # 7 FXI
    xₒ[8] = training_df[index, :XII]        # 8 FXII
    xₒ[9] = (1e-14)*SF                      # 9 FIIa
    model["initial_condition_vector"] = xₒ

    # what is the output array?
    Y = Array{Float64,1}(undef,5)
    Y[1] = training_df[index, :Lagtime]
    Y[2] = training_df[index, :Peak]
    Y[3] = training_df[index, Symbol("T.Peak")]
    Y[4] = training_df[index, :Max]
    Y[5] = training_df[index, :AUC]

    # setup parameters, and parameter bounds -
    if (isnothing(pₒ) == true)
        κ = [

            # default: hand fit set -
            0.061   0.01 10.0   ; # 1
            1.0     0.01 10.0   ; # 2
            1.0     0.01 10.0   ; # 3
            1.0     0.01 10.0   ; # 4
            1.0     0.01 10.0   ; # 5
            1.0     0.01 10.0   ; # 6
            1.0     0.01 10.0   ; # 7
            1.0     0.01 10.0   ; # 8
            0.70    0.01 10.0   ; # 9
            0.1     0.01 10.0   ; # 10
            0.045   0.01 10.0   ; # 11
            0.065   0.01 10.0   ; # 12
        ];

        # set default set as the start -
        pₒ = κ[:,1]
    end

    # what is my initial loss?
    lₒ = loss(pₒ, Y, model)

    for i ∈ 1:NTRIALS

        # generate a new candidate solution -
        Δᵢ = abs.(pₒ .* (1 .+ σ*rand(-1:1,12)))

        # compute the loss with this new solution -
        lᵢ = loss(Δᵢ, Y, model)

        # ok, is norm(lᵢ) < norm(lₒ)?
        if (norm(lᵢ) < norm(lₒ))

            # replace the "mom" -
            pₒ = Δᵢ

            # update the loss -
            lₒ = lᵢ

            @info "New soln for training set $(index) found at iteration $(i) with || lᵢ || = $(norm(lᵢ))"
        end
    end

    # run the sim w/the best parameters -
    # 1 - 9 : α vector
    model["α"] = pₒ[1:9]
    G = model["G"]
    idx = indexin(model, "FVIIa")   # 10
    G[idx, 4] = pₒ[10]
    idx = indexin(model, "AT")      # 11
    G[idx, 9] = pₒ[11]
    idx = indexin(model, "TFPI")    # 12
    G[idx, 1] = -1*pₒ[12]

    # run the model -
    (T,U) = evaluate(model)
    Xₘ = hcat(U...)
    Yₘ = model_output_vector(T,Xₘ[9,:]) # properties of the Thrombin curve

    return (pₒ, Yₘ, Y)
end

max_stars_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_issues_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_forks_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation


text:
# using Revise
using Test
using Logging
using SparseArrays
using Distributed
using LarSurf

@testset "Extract surface grid" begin
    segmentation = LarSurf.data234()
    verts, faces= LarSurf.get_surface_grid_per_voxel(segmentation; voxelsize_mm=[1,1,1])
    # verts are transposed, faces are in 2-D Array
    trifaces = LarSurf.triangulation(faces)

    # Visualization
    @test LarSurf.check_vf(verts, trifaces)

end


text:
function nlp_cvx_109_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - convex logarithmic objective
    # - binding nonlinear constraint
    # Variants
    #   010 - binding constraint (inflection point)
    #   011 - binding constraint (non-inflection point)

    model = Model(optimizer)

    @variable(model, x >= 0.00001, start = 0.1)
    @variable(model, y >= 0.00001, start = 0.1)

    @NLobjective(model, Max, log(x))
    @NLconstraint(model, (y - 2)^2 <= -x + 2)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, log(2), tol = objective_tol)
    return check_solution([x, y], [2, 2], tol = primal_tol)
end