text:
# auxiliary functions

using Distributions;
import Base.rand;

struct CategoricalSamplerNew <: Sampleable{Univariate,Discrete}
    mass::Vector{Float64}
    category::Vector{Float64}
end

function rand(s::CategoricalSamplerNew)
    s.category[rand(Categorical(s.mass))]
end

function genScenario_old(pDistr)
    # generate a disruption time
    tSupport = [i for i in keys(pDistr.tDistrn)];
    tProb = [pDistr.tDistrn[i] for i in tSupport];
    tDistrObj = Categorical(tProb);
    t = rand(tDistrObj);

    # generate a disruption location
    ωSupport = [i for i in keys(pDistr.ωDistrn)];
    ωProb = [pDistr.ωDistrn[i] for i in ωSupport];
    ωDistrObj = Categorical(ωProb);
    ω = rand(ωDistrObj);
    return tSupport[t],ωSupport[ω];
end

function genScenario(pDistr)
    # generate a disruption time
    tSupport = [i for i in keys(pDistr.tDistrn)];
    tProb = [pDistr.tDistrn[i] for i in tSupport];
    tDistrObj = Categorical(tProb);
    t = rand(tDistrObj);

    # generate a disruption location
    ωSupport = [i for i in keys(pDistr.ωDistrn)];
    ωProb = [pDistr.ωDistrn[i] for i in ωSupport];
    ωDistrObj = Categorical(ωProb);
    ω = rand(ωDistrObj);
    scenInd = ωSupport[ω];
    return tSupport[t],pDistr.ωDict[scenInd],pDistr.ωτ[scenInd];
end

function modifyOmega(pDistr,hardComp)
    ωDistrNew = Dict();
    ωDictNew = Dict();
    ωτNew = Dict();
    for iHard in hardComp
        releaseProb += pDistr.ωDistrn[iHard];
    end
    avgNo = length(values(pDistr.ωDistrn)) - length(hardComp);
    for i in keys(pDistr.ωDistrn)
        if !(i in hardComp)
            # if it is not the hardened component
            ωDistrNew[i] = pDistr.ωDistrn[i] + releaseProb/avgNo;
            ωDictNew[i] = pDistr.ωDict[i];
            ωτNew[i] = pDistr.ωτ[i];
        end
    end
    pDistrNew = probDistrn(pDistr.tDistrn,ωDistrNew,ωDictNew);
    return pDistrNew;
end

function modifyT(pDistr,λD,T)
    # modify the time distribution using the given λD
    tDistrnNew = Dict();
    for t in 1:(T-1)
        tDistrnNew[t] = exp(-λD*(t - 1)) - exp(-λD*t);
    end
    # probability of the time T+1
    tDistrnNew[T] = 1 - sum([tDistrnNew[t] for t in 1:(T-1)]);
    pDistrNew = probDistrn(tDistrnNew,pDistr.ωDistrn,pDistr.ωDict,pDistr.ωτ);
    return pDistrNew;
end

function calObj(mp, td, τ, T, fData, bData, pDistr, sp, lpp, lqp, lpm, lqm, θ)
    # function to formulate the objective function for f_D
    Ω = [ω for ω in keys(pDistr.ωDistrn)];

    objExpr = @expression(mp, 0);
    for tp in 1:maximum(keys(pDistr.tDistrn))
        dExpr = @expression(mp, 0);
        if tp <= T - (td + τ)
            for t in td:(tp + td + τ - 1)
                for i in fData.genIDList
                    # add generator cost
                    if fData.cp[i].n == 3
                        dExpr += fData.cp[i].params[1]*(sp[i,t]^2) + fData.cp[i].params[2]*sp[i,t];
                    elseif fData.cp[i].n == 2
                        dExpr += fData.cp[i].params[1]*sp[i,t];
                    end
                end
                # add load shed cost
                dExpr += fData.cz*(sum(lpp[i,t] + lqp[i,t] + lpm[i,t] + lqm[i,t] for i in fData.IDList));
            end
            objExpr += pDistr.tDistrn[tp]*(dExpr + sum(pDistr.ωDistrn[ω]*θ[tp + td + τ,ω] for ω in Ω));
        else
            for t in td:T
                for i in fData.genIDList
                    # add generator cost
                    if fData.cp[i].n == 3
                        dExpr += fData.cp[i].params[1]*(sp[i,t]^2) + fData.cp[i].params[2]*sp[i,t];
                    elseif fData.cp[i].n == 2
                        dExpr += fData.cp[i].params[1]*sp[i,t];
                    end
                end
                # add load shed cost
                dExpr += fData.cz*(sum(lpp[i,t] + lqp[i,t] + lpm[i,t] + lqm[i,t] for i in fData.IDList));
            end
            objExpr += pDistr.tDistrn[tp]*dExpr;
        end
    end
    return objExpr;
end

function calObjDet(mp, T, fData, bData, sp, lpp, lqp, lpm, lqm, u)
    # function to formulate the objective function for deterministic case

    objExpr = @expression(mp, sum(bData.cost[i]*u[i] for i in bData.IDList) +
        fData.cz*sum(sum(lpp[i,t] + lqp[i,t] + lpm[i,t] + lqm[i,t] for i in fData.IDList) for t in 1:T));
    for t in 1:T
        for i in fData.genIDList
            # add generator cost
            if fData.cp[i].n == 3
                objExpr += fData.cp[i].params[1]*(sp[i,t]^2) + fData.cp[i].params[2]*sp[i,t];
            elseif fData.cp[i].n == 2
                objExpr += fData.cp[i].params[1]*sp[i,t];
            end
        end
    end
    return objExpr;
end

function calObj1(mp, T, fData, bData, pDistr, sp, lpp, lqp, lpm, lqm, θ, u)
    # function to formulate the objective function for f_1
    Ω = [ω for ω in keys(pDistr.ωDistrn)];

    objExpr = @expression(mp, sum(bData.cost[i]*u[i] for i in bData.IDList));
    for tp in 1:maximum(keys(pDistr.tDistrn))
        dExpr = @expression(mp, 0);
        if tp < T
            for t in 1:tp
                for i in fData.genIDList
                    # add generator cost
                    if fData.cp[i].n == 3
                        dExpr += fData.cp[i].params[1]*(sp[i,t]^2) + fData.cp[i].params[2]*sp[i,t];
                    elseif fData.cp[i].n == 2
                        dExpr += fData.cp[i].params[1]*sp[i,t];
                    end
                end
                # add load shed cost
                dExpr += fData.cz*(sum(lpp[i,t] + lqp[i,t] + lpm[i,t] + lqm[i,t] for i in fData.IDList));
            end
            objExpr += pDistr.tDistrn[tp]*(dExpr + sum(pDistr.ωDistrn[ω]*θ[tp + 1,ω] for ω in Ω));
        else
            for t in 1:T
                for i in fData.genIDList
                    # add generator cost
                    if fData.cp[i].n == 3
                        dExpr += fData.cp[i].params[1]*(sp[i,t]^2) + fData.cp[i].params[2]*sp[i,t];
                    elseif fData.cp[i].n == 2
                        dExpr += fData.cp[i].params[1]*sp[i,t];
                    end
                end
                # add load shed cost
                dExpr += fData.cz*(sum(lpp[i,t] + lqp[i,t] + lpm[i,t] + lqm[i,t] for i in fData.IDList));
            end
            objExpr += pDistr.tDistrn[tp]*dExpr;
        end
    end
    return objExpr;
end

function calDualC1(T, fData, pDistr)
    # function to formulate the objective function for f_D
    Ω = [ω for ω in keys(pDistr.ωDistrn)];

    fDict = Dict();
    lDict = Dict();
    θDict = Dict();
    for tp in 1:maximum(keys(pDistr.tDistrn))
        if tp < T
            for t in 1:tp
                for i in fData.genIDList
                    if !((i,t) in keys(fDict))
                        fDict[i,t] = pDistr.tDistrn[tp];
                    else
                        fDict[i,t] += pDistr.tDistrn[tp];
                    end
                end
                for i in fData.IDList
                    if !((i,t) in keys(lDict))
                        lDict[i,t] = pDistr.tDistrn[tp]*fData.cz;
                    else
                        lDict[i,t] += pDistr.tDistrn[tp]*fData.cz;
                    end
                end
            end
            for ω in Ω
                θDict[tp + 1,ω] = pDistr.tDistrn[tp]*pDistr.ωDistrn[ω];
            end
        else
            for t in 1:T
                for i in fData.genIDList
                    if !((i,t) in keys(fDict))
                        fDict[i,t] = pDistr.tDistrn[tp];
                    else
                        fDict[i,t] += pDistr.tDistrn[tp];
                    end
                end
                for i in fData.IDList
                    if !((i,t) in keys(lDict))
                        lDict[i,t] = pDistr.tDistrn[tp]*fData.cz;
                    else
                        lDict[i,t] += pDistr.tDistrn[tp]*fData.cz;
                    end
                end
            end
        end
    end
    return fDict,lDict,θDict;
end

function calCostF(costn, currentSol, T, fData, nowT, disT)
    # function to calculate costs for forward pass
    costn += sum(sum(fData.cz*(abs(currentSol.lp[i,t]) + abs(currentSol.lq[i,t])) for i in fData.IDList) for t in nowT:(disT - 1));
    for t in nowT:(disT - 1)
        for i in fData.genIDList
            # add generator cost
            if fData.cp[i].n == 3
                costn += fData.cp[i].params[1]*(currentSol.sp[i,t]^2) + fData.cp[i].params[2]*currentSol.sp[i,t];
            elseif fData.cp[i].n == 2
                costn += fData.cp[i].params[1]*currentSol.sp[i,t];
            end
        end
    end
    return costn;
end

function calDualC(td, τ, T, fData, pDistr)
    # function to formulate the objective function for f_D
    Ω = [ω for ω in keys(pDistr.ωDistrn)];

    fDict = Dict();
    lDict = Dict();
    θDict = Dict();
    for tp in 1:maximum(keys(pDistr.tDistrn))
        if tp <= T - (td + τ)
            for t in td:(tp + td + τ - 1)
                for i in fData.genIDList
                    if !((i,t) in keys(fDict))
                        fDict[i,t] = pDistr.tDistrn[tp];
                    else
                        fDict[i,t] += pDistr.tDistrn[tp];
                    end
                end
                for i in fData.IDList
                    if !((i,t) in keys(lDict))
                        lDict[i,t] = pDistr.tDistrn[tp]*fData.cz;
                    else
                        lDict[i,t] += pDistr.tDistrn[tp]*fData.cz;
                    end
                end
            end
            for ω in Ω
                θDict[tp + td + τ,ω] = pDistr.tDistrn[tp]*pDistr.ωDistrn[ω];
            end
        else
            for t in td:T
                for i in fData.genIDList
                    if !((i,t) in keys(fDict))
                        fDict[i,t] = pDistr.tDistrn[tp];
                    else
                        fDict[i,t] += pDistr.tDistrn[tp];
                    end
                end
                for i in fData.IDList
                    if !((i,t) in keys(lDict))
                        lDict[i,t] = pDistr.tDistrn[tp]*fData.cz;
                    else
                        lDict[i,t] += pDistr.tDistrn[tp]*fData.cz;
                    end
                end
            end
        end
    end
    return fDict,lDict,θDict;
end

function simuPath(T,pDistr)
    pathList = [];
    nowT = 1;
    while nowT <= T
        tp,ωd,τω = genScenario(pDistr);
        push!(pathList, (tp,ωd,τω));
        if nowT == 1
            nowT += tp;
            nowT = min(nowT, T + 1);
        else
            nowT += tp + τω;
            nowT = min(nowT, T + 1);
        end
    end
    return pathList;
end

function pathSimu_cover(N,T,pDistr,genCutsCount,iterNo)
    # maximize the coverage from the sample
    pathSet = [];
    pathTimeList = [];
    for i in 1:(10*N)
        pathList = simuPath(T,pDistr);
        push!(pathSet,pathList);
        pathTimes = [];
        tNow = 1;
        for j in 1:length(pathList)
            if tNow == 1
                tNow += pathList[j][1];
            else
                tNow += pathList[j][3] + pathList[j][1];
            end
            if tNow < T
                push!(pathTimes, tNow);
            end
        end
        push!(pathTimeList, pathTimes);
    end

    sp = Model(solver = GurobiSolver(OutputFlag = 0));
    @variable(sp, x[i in 1:length(pathSet)], Bin);
    @constraint(sp, sum(x[i] for i in 1:length(pathSet)) == N);
    @objective(sp, Max, sum(sum((N*iterNo - genCutsCount[t])*x[i] for t in 2:T if t in [pathTimeList[i][j] for j in 1:length(pathTimeList[i])])
        for i in 1:length(pathSet)));
    solve(sp);
    pathSetsel = Dict();
    pathTimesel = Dict();
    iCount = 0;
    for i in 1:length(pathSet)
        if getvalue(sp[:x][i]) == 1
            iCount += 1;
            pathSetsel[iCount] = pathSet[i];
            pathTimesel[iCount] = pathTimeList[i];
        end
    end
    return pathSetsel,pathTimesel;
end

function pathSimu_cover_last(N,T,pDistr)
    # maximize the coverage for the terminal stage problem from the sample
    pathSet = [];
    pathTimeList = [];
    selNo = 0;
    while selNo < N
        pathList = simuPath(T,pDistr);
        pathTimes = [];
        tNow = 1;
        τList = false;
        for j in 1:length(pathList)
            if tNow == 1
                tNow += pathList[j][1];
            else
                tNow += pathList[j][3] + pathList[j][1];
            end
            if tNow < T
                push!(pathTimes, tNow);
                if tNow >= T - pathList[j][3]
                    τList = true;
                end
            end
        end
        if τList
            push!(pathSet,pathList);
            push!(pathTimeList, pathTimes);
            selNo += 1;
        end
    end
    return pathSet,pathTimeList;
end

# initialize the cutData in every core
function cutIni(cutData)
    global cutDict = cutData;
end

# update the cutData in every core with the new cuts
function cutUpdate(td,Ω,paraSet,cutCurrentData)
    # cutCurrentData is a list
    for ω in Ω
        itemInd = 0;
        for item in paraSet
            itemInd += 1;
            if item[1] == ω
                if (cutCurrentData[itemInd].solStatus == :Optimal)
                    if (td,ω) in keys(cutDict)
                        push!(cutDict[td,ω],cutCurrentData[itemInd]);
                    else
                        cutDict[td,ω] = [cutCurrentData[itemInd]];
                    end
                end
            end
        end
    end
end

function yzCal(y,bData,loc)
    noPiece = length(bData.ηα[loc]);
    fv = Inf;
    for n in 1:noPiece
        if y*bData.ηα[loc][n]+bData.ηβ[loc][n] < fv
            fv = y*bData.ηα[loc][n]+bData.ηβ[loc][n];
        end
    end
    return fv;
end

# change the cost of generation by a multiplier
function changeCost(fDatal, bDatal, cmulti)
    fDataLocal = deepcopy(fDatal);
    bDataLocal = deepcopy(bDatal);
    for i in fDataLocal.genIDList
        for nIter in 1:fDataLocal.cp[i].n
            fDataLocal.cp[i].params[nIter] = fDataLocal.cp[i].params[nIter]*cmulti^(fDataLocal.cp[i].n - nIter);
        end
    end
    fDataLocal.cz = fDataLocal.cz*cmulti;

    for i in bDataLocal.IDList
        bDataLocal.cost[i] = bDataLocal.cost[i]*cmulti;
    end

    global fData = fDataLocal;
    global bData = bDataLocal;
end

# break a component
function breakComponent(fDatal, bItem, bType)
    fDataLocal = deepcopy(fDatal);
    if bType == "g"
        # generators
        gListTemp = copy(fDataLocal.genIDList);
        deleteat!(gListTemp,findfirst(fDataLocal.genIDList,bItem));
        fDataLocal.genIDList = gListTemp;
    elseif bType == "l"
        # lines
        lineListTemp = deepcopy(fDataLocal.brList);
        for k in lineIDList
            if ((k[1],k[2]) == bItem) | ((k[2],k[1]) == bItem)
                deleteat!(lineListTemp,findfirst(fDataLocal.brList,k));
            end
        end
        fDataLocal.brList = lineListTemp;
    end
    global fData = fDataLocal;
end

# Add τ to the last element of pathDict
function addTau(pathDict,τ)
    newpathDict = Dict();
    for i in keys(pathDict)
        newPath = [];
        for item in pathDict[i]
            push!(newPath,(item[1],item[2],τ));
        end
        newpathDict[i] = newPath;
    end
    return newpathDict;
end

# get the index of scenarios
function reverseScen(pathDict,τ,pDistr)
    newpathDict = Dict();
    for i in keys(pathDict)
        newPath = [];
        for item in pathDict[i]
            # find the index of scenario
            jInd = 0;
            for j in keys(pDistr.ωDict)
                if item[2] in pDistr.ωDict[j]
                    jInd = j;
                end
            end
            push!(newPath,(item[1],pDistr.ωDict[jInd],pDistr.ωτ[jInd]));
        end
        newpathDict[i] = newPath;
    end
    return newpathDict;
end


text:
using Turing: Tracker

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

function turing_inference(
    prob::DiffEqBase.DEProblem,
    alg,
    t,
    data,
    priors;
    likelihood_dist_priors = [InverseGamma(2, 3)],
    likelihood = (u,p,t,σ) -> MvNormal(u, σ[1]*ones(length(u))),
    num_samples=1000, sampler = Turing.NUTS(0.65),
    syms = [Turing.@varname(theta[i]) for i in 1:length(priors)],
    sample_u0 = false,
    save_idxs = nothing,
    progress = false,
    kwargs...,
)
    N = length(priors)
    Turing.@model mf(x, ::Type{T} = Float64) where {T <: Real} = begin
        theta = Vector{T}(undef, length(priors))
        for i in 1:length(priors)
            theta[i] ~ NamedDist(priors[i], syms[i])
        end
        σ = Vector{T}(undef, length(likelihood_dist_priors))
        for i in 1:length(likelihood_dist_priors)
            σ[i] ~ likelihood_dist_priors[i]
        end
        nu = length(prob.u0)
        u0 = convert.(T, sample_u0 ? theta[1:nu] : prob.u0)
        p = convert.(T, sample_u0 ? theta[(nu + 1):end] : theta)
        _saveat = isnothing(t) ? Float64[] : t
        sol = concrete_solve(prob, alg, u0, p; saveat = _saveat, progress = progress, save_idxs = save_idxs, kwargs...)
        failure = size(sol, 2) < length(_saveat)

        if failure
            @logpdf() = T(0) * sum(x) + T(-Inf)
            return
        end
        if ndims(sol) == 1
            x ~ likelihood(sol[:], theta, Inf, σ)
        else
            for i = 1:length(t)
                x[:, i] ~ likelihood(sol[:, i], theta, sol.t[i], σ)
            end
        end
        return
    end

    # Instantiate a Model object.
    model = mf(data)
    chn = sample(model, sampler, num_samples; progress = progress)
    return chn
end

ext: jl


text:
#import Pkg; Pkg.add("Plots")
#import Pkg; Pkg.add("DataFrames")
#import Pkg; Pkg.add("CSV")

import Random

using Printf
using Statistics
using CSV, DataFrames
using DelimitedFiles

const ext = true
const parset = "figure1"
const seed = parse(Int64,ARGS[1])

if ext
    const mfpt_name = string("mfpt",parset)

else
    const mfpt_name = string("mfpt",parset,"_noext")
end

# note that switch_v values are from the mean-field model
# so forces and the switch_v must be scaled by 2

println("TEST0")

if mfpt_name == "mfpt1a"
	const p1=4
	const gamma=0.322
	const zeta = .2
	const nX = 100 # 12-100  on each side
	const nY = 100 # 12-100  on each side

	const alpha = 5 # rate constants 1/s
	const beta = 200

	const A = 5 # nm displacement of newly attached crossbridge
	const B = 5.05

	const T0 = 0#10/(alpha+beta) # start time
	const Tf = 50 # final time
	const V0 = 0 # nm/s shortening velocity for t > Tstart

	const switch_v = 61

elseif mfpt_name == "mfpt1b"
	const p1=4/5
	const gamma=0.322
	const zeta = .2
	const nX = 500 # 12-100  on each side
	const nY = 500 # 12-100  on each side

	const alpha = 5 # rate constants 1/s
	const beta = 200

	const A = 5 # nm displacement of newly attached crossbridge
	const B = 5.05

	const T0 = 0#10/(alpha+beta) # start time
	const Tf = 50 # final time
	const V0 = 0 # nm/s shortening velocity for t > Tstart

	const switch_v = 61

elseif mfpt_name == "mfpt1a_noext"

	const p1=4
	const gamma=0.322
	const zeta = .4
	const nX = 100 # 12-100  on each side
	const nY = 100 # 12-100  on each side

	const alpha = 10 # rate constants 1/s
	const beta = 150

	const A = 5 # nm displacement of newly attached crossbridge
	const B = 5.001

	const T0 = 0#10/(alpha+beta) # start time
	const Tf = 50 # final time
	const V0 = 0 # nm/s shortening velocity for t > Tstart

	const switch_v = 100

elseif mfpt_name == "mfpt1b_noext"

	const p1=4
	const gamma=0.322
	const zeta = .4
	const nX = 100 # 12-100  on each side
	const nY = 100 # 12-100  on each side

	const alpha = 10 # rate constants 1/s
	const beta = 150

	const A = 5 # nm displacement of newly attached crossbridge
	const B = 5.001

	const T0 = 0#10/(alpha+beta) # start time
	const Tf = 50 # final time
	const V0 = 0 # nm/s shortening velocity for t > Tstart

	const switch_v = 100

elseif mfpt_name == "mfpt2a_noext"

	const p1=4
	const gamma=0.322
	const zeta = .7
	const nX = 100 # 12-100  on each side
	const nY = 100 # 12-100  on each side

	const alpha = 14 # rate constants 1/s
	const beta = 300

	const A = 5 # nm displacement of newly attached crossbridge
	const B = 5.01

	const T0 = 0#10/(alpha+beta) # start time
	const Tf = 50 # final time
	const V0 = 0 # nm/s shortening velocity for t > Tstart

	const switch_v =  37

elseif mfpt_name == "mfpt2a"

	const p1=4
	const gamma=0.322
	const zeta = .7
	const nX = 100 # 12-100  on each side
	const nY = 100 # 12-100  on each side

	const alpha = 14 # rate constants 1/s
	const beta = 300

	const A = 5 # nm displacement of newly attached crossbridge
	const B = 5.01

	const T0 = 0#10/(alpha+beta) # start time
	const Tf = 50 # final time
	const V0 = 0 # nm/s shortening velocity for t > Tstart

	const switch_v =  36

elseif mfpt_name == "mfptfigure1"

	const p1=4
	const gamma=0.322
	const zeta = .4
	const nX = 100 # 12-100  on each side
	const nY = 100 # 12-100  on each side

	const alpha = 14 # rate constants 1/s
	const beta = 126

	const A = 5 # nm displacement of newly attached crossbridge
	const B = 5.05

	const T0 = 0#10/(alpha+beta) # start time
	const Tf = 25 # final time
	const V0 = 0 # nm/s shortening velocity for t > Tstart

	const switch_v =  121

end

println("p1=",p1," gamma=",gamma," zeta=",zeta," nX=",nX," nY=",nY," alpha=",alpha," beta=",beta," A=",A," B=",B," T0=",T0," Tf=",Tf," V0=",V0," switch_v=",switch_v," mfpt_name=",mfpt_name)

track_positions = false
attach_time = false
use_storage = false
use_last = false

if attach_time
    track_positions = true
end

function run_agents(dt::Float64)

    Random.seed!(seed)

    println("Initializing with seed=",seed)

    dt::Float64 = dt # time step

    a1 = zeros(UInt8,nX)#{Int64}(undef,nX) # binary variable for motors
    x1 = zeros(Float64,nX)#Array{Float64}(undef,nX) # position for each motor

    a2 = zeros(UInt8,nY)#Array{Int64}(undef,nY) # binary variable for motors
    x2 = zeros(Float64,nY)#Array{Float64}(undef,nY) # position for each motor

    #println(a1[change1])

    U1::Float64 = 0. # total attached
    U2::Float64 = 0. # total attached

    F1::Float64 = 0. # force
    F2::Float64 = 0. # force

    # position and velocity
    #V = V0

    #V = zeros(Float64,convert(Int64,floor(Tf/dt)+3))
    V::Float64 = V0

    switch_times = zeros(Float64,0)
    #println(switch_times)

    side::UInt8 = 0
    k::UInt64 = 1
    t::Float64 = 0.

    randsize::UInt16 = 50

    j::UInt64 = 1
    #ran1 = rand(randsize,nX)
    #ran2 = rand(randsize,nY)

	betadt::Float64 = beta*dt
	alphadt::Float64 = alpha*dt
	p1g::Float64 = p1*gamma

    while t < Tf

        t = (k-1)*dt

        #if !ext
        #    dx1 = 0
        #    dx2 = 0
        #else
        #    dx1 = V*dt
        #    dx2 = -V*dt
        #end

		if V >= 0
			dx1 = V*dt
			dx2 = -ext*V*dt
		else
			dx1 = ext*V*dt
			dx2 = -V*dt
		end

        for i = 1:length(a1)
            if a1[i]>0

                x1[i] = x1[i] + dx1

               if ((rand() < betadt) | (x1[i] > B))
               #if ((ran1[j,i] < (beta*dt)) | (x1[i] > B))
                    a1[i] = 0
                    x1[i] = 0
                end
            else
                if (rand() < alphadt)# & (V < 0)
                #if (ran1[j,i] < (alpha*dt))
                    a1[i] = 1
                    x1[i] = A
                end
            end
        end

        for i = 1:length(a2)
            if a2[i] > 0

                x2[i] = x2[i] + dx2

                if ((rand() < betadt) | (x2[i] > B))
                #if ((ran2[j,i] < (beta*dt)) | (x2[i] > B))
                    a2[i] = 0
                    x2[i] = 0
                end
            else
                if (rand() < alphadt)# & (V > 0)
                #if (ran2[j,i] < (alpha*dt))
                    a2[i] = 1
                    x2[i] = A
                end
            end
        end

        # total motor number at time step j
		U1 = 0
		U2 = 0

		for i = 1:length(a1)
			U1 = U1 + a1[i]
		end

		for i = 1:length(a2)
			U2 = U2 + a2[i]
		end
        #U1 = sum(a1)
        #U2 = sum(a2)

        # net force 1
        # x*p1*gamma

		F1 = 0
		F2 = 0
		for i = 1:length(x1)
			F1 = F1 + x1[i]
		end

		for i = 1:length(x2)
			F2 = F2 + x2[i]
		end

		F1 = F1*p1g
		F2 = F2*p1g

        #F1 = sum(x1)*p1g
        #F2 = sum(x2)*p1g

        #V += dt*(-F1+F2-zeta*V)/eps
        #V[k+1] = V[k] + dt*(F2-F1-zeta*V[k])/eps
		V = (-F1+F2)/zeta

        if (side == 0) & (V >= switch_v)
        #if (side == 0) & (V[k+1] >= switch_v)
            side = 1
            append!(switch_times,t)
        end

        if (side == 1) & (V <= -switch_v)
        #if (side == 1) & (V[k+1] <= -switch_v)
            side = 0
            append!(switch_times,t)
        end

        #print(V)
        j += 1
        k += 1
    end

    return V,switch_times

end

#using ProfileView

#@profview run_agents(0.01)
#@profview run_agents(0.0001)

#for dt in [0.001 0.0001 0.00001 0.000001 0.0000001]
for dt in [0.000003]
#for dt in [0.0000075 0.0000005 0.00000025 0.0000001]
#for dt in [0.00000005 0.00000002]

	mfpt_dir = string(mfpt_name)

	if ~isdir(mfpt_dir)
		#println("WARNING: DIRECTORY",mfpt_dir," DOES NOT EXIST. PLEASE CREATE.")
		println("creating directory ",mfpt_dir)
		mkdir(mfpt_dir)
	end

	V,out = run_agents(dt)
	times = diff(out[out.>5])

	if false
		plt = plot!(V,show=true)
		gui(plt)
		savefig(plt,"temp.png")
	end

	#fname_tint = "mfpt/TEMP_time_intervals_Tf=%d_nX=%d_nY=%d_seed=%d_dt=%g.txt"%(Tf,nX,nY,seed,round_to_n(dt))

	#@sprintf("%.g",mfpt_dir)

	open(@sprintf("%s/time_intervals_Tf=%d_nX=%d_nY=%d_seed=%d_dt=%.1e_alpha=%g_beta=%g.txt",mfpt_dir,Tf,nX,nY,seed,dt,alpha,beta),"w") do io
		writedlm(io,times)
	end

	#println(times)
	println("__MFPT__=",mean(times),", total switches=",length(times),", dt=",dt)

end

#end
#CSV.write(fname,diff(out))

max_stars_repo_name: youngmp/NoisyMotors

max_issues_repo_name: youngmp/NoisyMotors

max_forks_repo_name: youngmp/NoisyMotors


text:
model,obj = Dynamics.dubinscar
solver = Solver(model,obj,N=21)
n,m,N = get_sizes(solver)
obj.cost.Q .= Diagonal(I,n)
obj.cost.R .= Diagonal(I,m)*2
obj.cost.H .= ones(m,n)

X = rand(n,N)
U = rand(m,N-1)
dt = solver.dt
dx = zeros(n,n)
du = zeros(n,m)
x = X[:,1]
u = U[:,1]
solver.Fd(dx,du,x,u)
dz = zeros(n,n+m)
z = [x;u]
solver.Fd(dz,z)
dz == [dx du]

ind_z = create_partition([n,m],(:x,:u))
ind2_z = create_partition2([n,m],(:x,:u))
f = model.f
fd! = rk4(f,solver.dt)
fd_aug!(xdot,z) = fd!(xdot,view(z,ind_z.x),view(z,ind_z.u))
f_aug!(xdot,z) = f(xdot,view(z,ind_z.x),view(z,ind_z.u))
fd_aug2!(dS,S::Array) = begin
    k1 = k2 = k3 = k4 = zero(S)
    f_aug!(k1,S);         k1 *= dt;
    f_aug!(k2,S + k1/2); k2 *= dt;
    f_aug!(k3,S + k2/2); k3 *= dt;
    f_aug!(k4,S + k3);    k4 *= dt;
    copyto!(dS, (S + (k1 + 2*k2 + 2*k3 + k4)/6)[1:n])
end
xdot = zeros(n)
fd_aug!(xdot,z)
xdot2 = zeros(n)
fd_aug2!(xdot2,z)
xdot2 == xdot

Fd!(dz,z) = ForwardDiff.jacobian!(dz,fd_aug!,xdot,z)
Fd2!(dz,z) = ForwardDiff.jacobian!(dz,fd_aug2!,xdot,z)
dz2 = zero(dz)
Fd!(dz,z)
Fd2!(dz2,z)
dz == dz2 == [dx du]
@btime Fd!($dz,$z)
@btime Fd2!($dz,$z)
@btime solver.Fd($dz,$z)

function Fd_split!(dx,du,x,u)
    z[ind_z.x] = x
    z[ind_z.u] = u
    Fd!(dz,z)
    copyto!(dx,dz[ind2_z.xx])
    copyto!(du,dz[ind2_z.xu])
end
dx2 = zero(dx)
du2 = zero(du)
solver.Fd(dx,du,x,u)
Fd_split!(dx2,du2,x,u)
dx == dx2
du == du2
@btime solver.Fd($dx,$du,$x,$u)
@btime Fd_split!($dx,$du,$x,$u)


text:
using .Pluto

@info "You can use the notebooks in `TropicalTensors` now by typing, e.g.
  ∘ `TropicalTensors.notebook(\"spinglass\")`, solving square lattice spinglass using a quantum simulator.
  ∘ `TropicalTensors.notebook(\"ising_and_2sat\")`, solving Ising spinglass and 2-SAT counting using tensor network contraction.
"

"""
    notebook(which; dev=false, kwargs...)

Open a notebook, the first argument can be
* "spinglass": solving spinglass model with Yao.
* "randomgraph": solving randomgraph model with tensor network contraction.
"""
function notebook(which; dev=false, kwargs...)
    src = project_relative_path("notebooks", "$which.jl")
    if dev
        dest = src
    else
        dest = tempname()
        cp(src, dest)
        chmod(dest, 0o664)
    end
    Pluto.run(; notebook=dest, project=project_relative_path(), kwargs...)
end


text:
function bayesreg(y,x)
  n = length(y)
  k = size(x,2)
  X = [ones(n) x]
  # uninformative analytical
  iXX = inv(X'*X)
  bhat = iXX*X'*y  # bhat = X \ y
  s2hat = (y - X*bhat)'*(y - X*bhat)/n
  Vb = s2hat.*iXX
  seb = sqrt.(diag(Vb))

  println("linear model results:")
  println("coeff estimates ", bhat)
  println("standard errors ",seb)
  println("s^2 (eqn. variance) = ",s2hat)

  # compute R^2
  tss = sum((y .- mean(y)).^2)
  R2 = 1.0 - s2hat*n/tss[1]
  println("Rsquared = ",R2)
  yhat = X*bhat;

  return bhat, seb, s2hat, Vb, R2, yhat
end

# conjugate NIG prior regression:
"""
  bayesregNIG(y,x)

conjugate NIG prior linear regression function: Compute coeff. estimates, s.e's, equation σ, Rsquared

usage:

    bhat, seb, s2hat, Vb, R2, yhat = bayesregNIG(y,x; b0=0.,B0=100000.,a0=0.00001,d0=0.00001)
"""
function bayesregNIG(y,x; b0=0.0,B0=100000.0,a0=0.00001,d0=0.00001)
  n = length(y)
  X = [ones(n) x]
  k = size(X,2)
  if length(b0) < 2
    b0 = b0.*ones(2)
  end
  k = size(X,2)
  if size(B0,1) < 2
    B0 = B0.*Matrix(1.0I, m, m)
  end

  # informative analytical- Greenberg, ch. 4, (4.5), (4.6)
  iB0 = inv(B0) # B0\1  # better way to invert computationally
  B1 = inv(X'*X + iB0)
  b1 = B1*(X'*y + iB0*b0)  # bhat = X \ y
  a1 = a0 + n - k
  postssr = (y - X*b1)'*(y - X*b1)
  d1 = d0 + postssr
  s2hat = d1/a1
  Vb = s2hat.*B1
  seb = sqrt.(diag(Vb))

  println("linear model posterior results:")
  println("coeff estimates ", b1)
  println("standard errors ",seb)
  println("s^2 (eqn. variance) = ",s2hat)

  # compute R^2
  tss = sum((y .- mean(y)).^2)
  R2 = 1.0 - s2hat*n/tss
  println("Rsquared = ",R2)
  yhat = X*b1

  return b1, seb, s2hat, Vb, R2, yhat
end

"""
  blinreg(y,x)

linear regression function: Compute coeff. estimates, s.e's, odds, equation σ, R-squared

usage:

  b, seb, odds, s, R2 = linreg(y,x)
"""
function blinreg(y,x)
  # add intercept
  n = length(y)
  X = [ones(n) x]
  b = (X'*X) \ X'*y
  resids = y - X*b
  s2 = sum(resids.^2)/n
  covb = s2.*inv(X'X)
  seb = sqrt.(diag(covb))
  k = length(b)
  odds = zeros(k)
  for i in 1:k
    odds[i] = todds(b[i],seb[i],(n-k))
  end
  println("coeffs = ",round.(b,3))
  println(" s.e's = ",round.(seb,3))
  println(" odds = ",round.(odds,3))
  println("s^2 (eqn. variance) = ",round(s2,6))
  # compute R^2
  tss = sum((y .- mean(y)).^2)
  R2 = 1.0 - s2*n/tss
  println("Rsquared = ",round(R2,3))
  s = sqrt(s2)
  return b, seb, odds, s, R2;
end


text:
using Test, ApproxFun, SingularIntegralEquations
import ApproxFunBase: ∞, testbandedoperator, testfunctional, testblockbandedoperator, testraggedbelowoperator
import ApproxFunOrthogonalPolynomials: JacobiZ
import SingularIntegralEquations: testsies, testsieeval, stieltjesmoment, Directed, _₂F₁, ⁺, ⁻

@testset "Log Kernel" begin
    @testset "Segment" begin
        testsieeval(Jacobi(0,0))

        a=1.0;b=2.0+im
        d=Segment(a,b)
        z=Fun(d)
        f=real(exp(z)/(sqrt(z-a)*sqrt(b-z)))
        S=space(f)
        testsies(S)

        a=1.0;b=2.0+im
        d=Segment(a,b)
        z=Fun(d)
        f=real(exp(z)*(sqrt(z-a)*sqrt(b-z)))
        S=space(f)
        testsies(S)

        a=1.0;b=2.0+im
        d=Segment(a,b)
        z=Fun(d)
        f=real(exp(z)/(sqrt(z-a)*sqrt(b-z)))
        S=JacobiWeight(-0.5,-0.5,ChebyshevDirichlet{1,1}(d))
        # TODO: move to testsies
        H=OffSingularIntegral(S,Chebyshev(3..4),0)
        @test (H*f)(3.5) ≈ logkernel(f,3.5)

        H=OffSingularIntegral(S,Chebyshev(Segment(3,4.0+im)),0)
        @test (H*f)(3.5+0.5im) ≈ logkernel(f,3.5+0.5im)
    end

    @testset "Circle" begin
        d=Circle(0.2,3.0)
        S=Fourier(d)
        ζ=Fun(d)
        f=real(ζ+1/(ζ-0.1))
        #z=0.1+0.1im;@test linesum(logabs(ζ-z)*f) ≈ logkernel(f,z)*π
        #z=5.0+0.1im;@test linesum(logabs(ζ-z)*f) ≈ logkernel(f,z)*π

        d=Circle(0.2,3.0)
        S=Fourier(d)
        H=Hilbert(S,0)
        ζ=Fun(d)
        f=real(ζ+1/(ζ-0.1))
        z=0.2+3im;@test (H*f)(z) ≈ logkernel(f,z)
    end

    @testset "Legendre and Jacobi on interval" begin
        x=Fun()
        f=exp(x)

        @test isa(logkernel(f,2.0+im),Real)
        @test logkernel(f,2.0+im) ≈ sum(f*logabs(2.0+im-x))/π

        f=sqrt(1-x^2)*exp(x)
        @test logkernel(f,2.0+im) ≈ sum(f*logabs(2.0+im-x))/π

        x=Fun()
        f=(1-x)^0.1
        sp=space(f)
        @test logjacobimoment(sp.α,sp.β,2.0) ≈ sum((1-x)^sp.α*(1+x)^sp.β*logabs(2.0-x))

        @test logkernel(f,2.0) ≈ sum(f*logabs(2.0-x)/π)

        f=(1-x)^0.1*exp(x)

        @test stieltjes(f,2.0) ≈ sum(f/(2.0-x))

        @test logkernel(f,2.0) ≈ sum(f*logabs(2.0-x)/π)

        @test isa(logkernel(f,2.0+im),Real)
        @test logkernel(f,2.0+im) ≈ sum(f*logabs(2.0+im-x)/π)

        f=(1-x^2)^0.1*exp(x)
        sp=space(f)
        @test logjacobimoment(sp.α,sp.β,2.0) ≈ sum((1-x^2)^0.1*logabs(2.0-x))
        @test logkernel(f,2.0+im) ≈ sum(f*logabs(2.0+im-x)/π)

        f=(1-x)^(-0.1)*(1+x)^(-0.2)*exp(x)
        @test logkernel(f,2.0+im) ≈ sum(f*logabs(2.0+im-x)/π)
        @test isa(logkernel(f,2.0+im),Real)
    end

    @testset "other intervals" begin
        x = Fun()
        f = exp(x)
        z = 2+im
        @test logkernel(f,z) ≈ linesum(f*logabs(x-z))/π
        @test logkernel(f, -10.0) ≈ linesum(f*logabs(x+10))/π

        x = Fun()
        f = real(sqrt(1-x)*exp(x))
        z = 2+im
        @test logkernel(f,z) ≈ linesum(f*logabs(x-z))/π
        @test logkernel(f, -10.0) ≈ linesum(f*logabs(x+10))/π
        @test_broken logkernel(f, 0.1) ≈ logkernel(f, 0.1+eps()im)

        x = Fun(0.1..1)
        f = real(sqrt(1-x)*exp(x))
        z = 2+im
        @test logkernel(f,z) ≈ linesum(f*logabs(x-z))/π
        @test logkernel(f, -10.0) ≈ linesum(f*logabs(x+10))/π
        @test_broken logkernel(f, 0.2) ≈ logkernel(f, 0.2+eps()im)
    end

    @testset "Arc" begin
        a=Arc(0.,1.,0.,π/2)
        testsieeval(Legendre(a);posdirection=(-1-im))
    end

    @testset "Chebyshev singularities" begin
        for d in (ChebyshevInterval(), Interval(0,1), Segment(1+im, 2+3im))
            for k = 1:10
                u = Fun(JacobiWeight.(-0.5,-0.5,Ref(ChebyshevDirichlet{1,1}(d))), [zeros(4);1])
                x = Fun(d)
                for z in (20.0+im, 3.0)
                    ex = linesum(u*logabs(z-x))/π
                    @test logkernel(u, z) ≈ ex
                end
            end
        end

        for d in (ChebyshevInterval(), Interval(0,1))
            for k = 1:10
                u = Fun(JacobiWeight.(-0.5,-0.5,Ref(ChebyshevDirichlet{1,1}(d))), [zeros(4);1])
                x = Fun(d)
                for z in (20.0+im, 3.0)
                    ex_s = sum(u*log(z-x))
                    @test stieltjesintegral(u, z) ≈ ex_s
                end
            end
        end
    end
end


text:
using ElectronicStructure
using ElectronicStructure.LagrangeFunction

include("init_pot_specials.jl")

function test_main( ; method = "Emin_cg" )
  # LF parameters
  NN = [25, 25, 25]
  AA = [0.0, 0.0, 0.0]
  BB = [6.0, 6.0, 6.0]

  Npoints = prod(NN)

  # Initialize LF
  LF = init_LF3d_p( NN, AA, BB )
  ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h

  # Initialize G-vectors
  L = BB - AA
  Gv = GvectorsT( NN, diagm(L) )

  # Parameter for potential
  center = AA + 0.5*(BB-AA)
  # Potential
  ω = 2.0
  V_ionic = init_pot_harm_3d( LF, ω, center )

  Ncols = 4
  Focc = 2.0*ones(Ncols)

  if method == "Emin_cg_sparse"
    #
    #∇2 = get_Laplacian3d_kron(LF)
    ∇2 = build_nabla2(LF)
    precH = prec_mkl_ilu0( -0.5*∇2 + spdiagm(V_ionic) )
    #precH = speye(Npoints)  # for testing unpreconditioned code
    Energies, evecs, Potentials = KS_solve_Emin_pcg( LF, Gv, ∇2, precH,
                                    V_ionic, Focc, Ncols, verbose=true )
    evals = calc_evals( LF, ∇2, Potentials, evecs )
    #
  elseif method == "SCF"
    ∇2 = get_Laplacian3d_kron(LF)
    precH = prec_mkl_ilu0( -0.5*∇2 + spdiagm(V_ionic) )
    #precH = speye(Npoints)
    Energies, evals, evecs, Potentials =
          KS_solve_scf( LF, Gv, ∇2, precH, V_ionic, Focc, Ncols, verbose=true )
    #
  else
    #
    Energies, evecs, Potentials = KS_solve_Emin_cg( LF, Gv, V_ionic, Focc, Ncols, verbose=true )
    evals = calc_evals( LF, Potentials, evecs )
  end

  @printf("\nEigenvalues:\n")
  for i = 1:Ncols
    @printf("%8d %f\n", i, evals[i])
  end

  rho = calc_rho( Focc, evecs )
  @printf( "\nintegRho = %lf\n", sum(rho)*ΔV )

end

@code_native test_main()
@time test_main(method="Emin_cg_sparse")
#@time test_main(method="Emin_cg")
#@time test_main(method="SCF")


text:
function solve_modelARCDiagAbs(H, g, nlp_stop, PData::PDataFact, α::T) where {T}
    #printsyled("on est dans solve_modelARCDiagAbs \n", color = :cyan)
    # Solve the ARC subproblem once diagonalized into Δ using the norm |Δ|

    # Setup the subproblem to recover identity regularizing matrix
    # Should this be in pre-process so that PData.λ is indeed used as updated?
    # No No! Use PData.success to initialize PData.λ when true, and use it otherwise
    #
    # TO TEST to confirm it avoids ascent directions
    #
    ḡ = PData.g̃
    #printsyled("on a ḡ = $ḡ\n", color = :yellow)
    n_g = norm(ḡ)
    #printsyled("on a n_g = $n_g\n", color = :yellow)
    # ϵ =  1.0e-10
    ϵ = sqrt(eps(T)) * 100.0
    #printsyled("on a ϵ = $ϵ  et PData.Δ = $(PData.Δ)\n", color = :yellow)
    Γ2 = max(maximum(abs.(PData.Δ)), ϵ)
    #printsyled("on a Γ2 = $Γ2 \n", color = :yellow)
    Γ = sqrt(Γ2)
    #printsyled("on a Γ = $Γ \n", color = :yellow)
    Δ̄ = Γ .\ PData.Δ ./ Γ
    #printsyled("on a Δ̄  = $Δ̄  \n", color = :yellow)
    Δ = PData.Δ
    #printsyled("on a Δ = $Δ \n", color = :yellow)
    λmin = 0.0 #  Test the following
    #printsyled("λmin = $λmin \n", color = :yellow)
    if PData.success # ensure
        #printsyled("PData.success = $(PData.success) \n", color = :yellow)
        l_m, = findmin(Δ̄)
        #printsyled("l_m = $l_m \n", color = :yellow)
        λ = max(-l_m, 0.0)
        #printsyled("λ = $λ \n", color = :yellow)
        λmin = max(ϵ, λ + ϵ * (1.0 + λ))
        #printsyled("λmin = $λmin \n", color = :yellow)
    else
        #printsyled("PData.success = $(PData.success) \n", color = :yellow)
        λmin = PData.λ
    end
    d̄ = -(Δ .+ λmin * Γ2) .\ ḡ
    #printsyled("d̄ = $d̄ \n", color = :yellow)
    seuil_bar = norm(d̄) / α

    #printsyled("on est a seuil_bar = $seuil_bar \n", color = :yellow)
    # Solve the subproblem (Δ̄ + λ I) d̄ = -ḡ such that λ = ||d̄||/α
    d̄, λ = solve_diag(λmin, Δ, ḡ, seuil_bar, α, ϵ, M = Γ2)
    # Transform back d_s into d
    d̃ = d̄
    #d̃ = Γ .\ d̄
    d = AInv(PData, d̃)

    #try assert((PData.g̃ + 0.5*PData.Δ .* d̃)⋅d̃ <= 0.0)  catch  @bp  end
    #println("*******SolveModelDiagAbs:  PData.g̃⋅d̃ = $(PData.g̃⋅d̃), 0.5 d̃'PData.Δd̃ = $(0.5*(PData.Δ .* d̃)⋅d̃)")

    PData.λ = λ

    return d, λ
end


text:
import Dates
using Mocking

utc = FixedTimeZone("UTC")
warsaw = first(compile("Europe/Warsaw", tzdata["europe"]))
apia = first(compile("Pacific/Apia", tzdata["australasia"]))
midway = first(compile("Pacific/Midway", tzdata["australasia"]))

@testset "Construct ZonedDateTime / DateTime" begin
    # Constructing a ZonedDateTime from a DateTime and the reverse
    dt = DateTime(2019, 4, 11, 0)
    zdt = ZonedDateTime(dt, warsaw)
    @test DateTime(zdt) == DateTime(2019, 4, 11, 0)
    @test DateTime(zdt, UTC) == DateTime(2019, 4, 10, 22)

    # Converting between ZonedDateTime and DateTime isn't possible as it isn't lossless.
    @test_throws MethodError convert(DateTime, zdt)
    @test_throws MethodError convert(ZonedDateTime, dt)
end

@testset "Construct Date / ZonedDateTime" begin
    date = Date(2018, 6, 14)
    zdt = ZonedDateTime(date, warsaw)
    @test Date(zdt) == Date(2018, 6, 14)
    @test Date(zdt, UTC) == Date(2018, 6, 13)

    # Converting between ZonedDateTime and Date isn't possible as it isn't lossless.
    @test_throws MethodError convert(Date, zdt)
    @test_throws MethodError convert(ZonedDateTime, date)
end

@testset "Construct Time" begin
    zdt = ZonedDateTime(2017, 8, 21, 0, 12, warsaw)
    @test Time(zdt) == Time(0, 12)
    @test Time(zdt, UTC) == Time(22, 12)

    # Converting between ZonedDateTime and Time isn't possible as it isn't lossless.
    @test_throws MethodError convert(Time, zdt)
end

@testset "now" begin
    dt = now(Dates.UTC)::DateTime
    zdt = now(warsaw)
    @test zdt.timezone == warsaw
    @test Dates.datetime2unix(DateTime(zdt, UTC)) ≈ Dates.datetime2unix(dt)
end

@testset "today" begin
    @test abs(today() - today(warsaw)) <= Dates.Day(1)
    @test today(apia) - today(midway) == Dates.Day(1)
end

@testset "todayat" begin
    @testset "current time" begin
        local zdt = now(warsaw)
        now_time = Time(zdt)
        @test todayat(now_time, warsaw) == zdt
    end

    @testset "ambiguous" begin
        local patch = @patch today(tz::TimeZone) = Date(1916, 10, 1)
        apply(patch) do
            @test_throws AmbiguousTimeError todayat(Time(0), warsaw)
            @test todayat(Time(0), warsaw, 1) == ZonedDateTime(1916, 10, 1, 0, warsaw, 1)
            @test todayat(Time(0), warsaw, 2) == ZonedDateTime(1916, 10, 1, 0, warsaw, 2)
        end
    end
end

# Changing time zones
dt = DateTime(2015, 1, 1, 0)
zdt_utc = ZonedDateTime(dt, utc; from_utc=true)
zdt_warsaw = ZonedDateTime(dt, warsaw; from_utc=true)

# Identical since ZonedDateTime is immutable
@test astimezone(zdt_utc, warsaw) === zdt_warsaw
@test astimezone(zdt_warsaw, utc) === zdt_utc

# ZonedDateTime to Unix timestamp (and vice versa)
@test TimeZones.zdt2unix(ZonedDateTime(1970, utc)) == 0
@test TimeZones.unix2zdt(0) == ZonedDateTime(1970, utc)

for dt in (DateTime(2013, 2, 13), DateTime(2016, 8, 11))
    local dt
    local zdt = ZonedDateTime(dt, warsaw)
    offset = TimeZones.value(zdt.zone.offset)   # Total offset in seconds
    @test TimeZones.zdt2unix(zdt) == datetime2unix(dt) - offset
end

@test isa(TimeZones.zdt2unix(ZonedDateTime(1970, utc)), Float64)
@test isa(TimeZones.zdt2unix(Float32, ZonedDateTime(1970, utc)), Float32)
@test isa(TimeZones.zdt2unix(Int64, ZonedDateTime(1970, utc)), Int64)
@test isa(TimeZones.zdt2unix(Int32, ZonedDateTime(1970, utc)), Int32)

@test TimeZones.zdt2unix(ZonedDateTime(1970, 1, 1, 0, 0, 0, 750, utc)) == 0.75
@test TimeZones.zdt2unix(Float32, ZonedDateTime(1970, 1, 1, 0, 0, 0, 750, utc)) == 0.75
@test TimeZones.zdt2unix(Int64, ZonedDateTime(1970, 1, 1, 0, 0, 0, 750, utc)) == 0
@test TimeZones.zdt2unix(Int32, ZonedDateTime(1970, 1, 1, 0, 0, 0, 750, utc)) == 0

# round-trip
zdt = ZonedDateTime(2010, 1, 2, 3, 4, 5, 999, utc)
round_trip = TimeZones.unix2zdt(TimeZones.zdt2unix(zdt))
@test round_trip == zdt

# millisecond loss
zdt = ZonedDateTime(2010, 1, 2, 3, 4, 5, 999, utc)
round_trip = TimeZones.unix2zdt(TimeZones.zdt2unix(Int64, zdt))
@test round_trip != zdt
@test round_trip == floor(zdt, Dates.Second(1))

# timezone loss
zdt = ZonedDateTime(2010, 1, 2, 3, 4, 5, warsaw)
round_trip = TimeZones.unix2zdt(TimeZones.zdt2unix(Int64, zdt))
@test round_trip == zdt
@test timezone(round_trip) != timezone(zdt)

# Julian dates
jd = 2457241.855
jd_zdt = ZonedDateTime(Dates.julian2datetime(jd), warsaw, from_utc=true)
@test TimeZones.zdt2julian(jd_zdt) == jd
@test TimeZones.zdt2julian(Int, jd_zdt) === floor(Int, jd)
@test TimeZones.zdt2julian(Float64, jd_zdt) === jd
@test TimeZones.julian2zdt(jd) == jd_zdt


text:
# Implementation is taken from [Flux.jl](https://github.com/FluxML/Flux.jl)
"""
    Chain(layers...)
Chain multiple layers / functions together, so that they are called in sequence
on a given input.
```julia
m = Chain(x -> x^2, x -> x+1)
m(5) == 26
m = Chain(Dense(input=10, output=5), Dense(input=5, output=2))
x = rand(10)
m(x) == m[2](m[1](x))
```
`Chain` also supports indexing and slicing, e.g. `m[2]` or `m[1:end-1]`.
`m[1:3](x)` will calculate the output of the first three layers.
"""
struct Chain{T<:Tuple}
  layers::T
   Chain(xs...) = new{typeof(xs)}(xs)
end

children(c::Chain) = c.layers
mapchildren(f, c::Chain) = Chain(f.(c.layers)...)

applychain(::Tuple{}, x) = x
applychain(fs::Tuple, x) = applychain(Base.tail(fs), first(fs)(x))

(c::Chain)(x) = applychain(c.layers, x)

Base.getindex(c::Chain, i::AbstractArray) = Chain(c.layers[i]...)
Base.getindex(c::Chain, i::Integer) = c.layers[i]
Base.getindex(c::Chain, ::Colon) = c


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ a450c12c-753d-11eb-06fc-29c5d907cb06
begin
	using Pkg
	Pkg.activate("MLJ_env", shared=true)

	# using MLJ
	using CSV
	using DataFrames
	using PlutoUI
	using Random
	using Test
	using Printf
	# using Plots
end

# ╔═╡ a4ce9686-7549-11eb-2172-2d47391d1260
include("./utils.jl")

# ╔═╡ 5b0c59ae-753d-11eb-1a8a-6fa0bf99a357
md"""
## Week 5: Lasso (Coordinate Descent)

In this notebook, you will implement your very own LASSO solver via coordinate descent. We will:

  - Write a function to normalize features
  - Implement coordinate descent for Lasso
  - Explore effects of L1 penalty

"""

# ╔═╡ 7e12df34-7549-11eb-30a5-b3099a8f77f6
md"""
### Load in house sales data
"""

# ╔═╡ 9e5e1060-7549-11eb-0418-55bcd0b8c757
sales = CSV.File("../../ML_UW_Spec/C02/data/kc_house_test_data.csv";
	header=true) |> DataFrame;

# ╔═╡ 2fa16676-754a-11eb-00cc-896d0e28b246
md"""
### Normalize features

In the house dataset, features vary wildly in their relative magnitude: sqft_living is very large overall compared to bedrooms, for instance. As a result, weight for sqft_living would be much smaller than weight for bedrooms. This is problematic because "small" weights are dropped first as l1_penalty goes up.

To give equal considerations for all features, we need to normalize features as discussed in the lectures: we divide each feature by its 2-norm so that the transformed feature has norm 1.
"""

# ╔═╡ 4e5d02be-754a-11eb-061f-41ab2f747c54
md"""
moved to `utils.jl`

```julia
using LinearAlgebra

function normalize_features(X::Matrix)
  n = size(X)[2]
  norms = zeros(eltype(X), n)'
  for ix ∈ 1:n
    norms[ix] = norm(X[:, ix])
  end
  (X ./ norms, norms)
end
```
"""

# ╔═╡ dc99317e-754a-11eb-0a8e-61ec5dff0507
begin
	M1 = [3. 6. 9.; 4. 8. 12.]
	features1, norms1 = normalize_features(M1)

	@test features1 == [0.6 0.6 0.6; 0.8 0.8 0.8]
	@test norms1 == [5.  10.  15.]  ## or [5.,  10.,  15.]
end

# ╔═╡ e9128af8-7555-11eb-1f4d-d1e17d70f4b3
typeof(features1), typeof(norms1)

# ╔═╡ 0d7f1fa0-754c-11eb-330b-d31a86f5566c
begin
	# 3 rows x 2 cols
	M2 = [3. 4.; 6. 5.; 9. 10.]
	features2, norms2 = normalize_features(M2)

	@test features2 ≈ [0.2672612419124244 0.3368607684266076;
		0.5345224838248488 0.42107596053325946;
		0.8017837257372732 0.8421519210665189]
	@test norms2 ≈ [11.224972160321824 11.874342087037917]
end

# ╔═╡ 984c778a-7557-11eb-21f8-dd2502fe7b73
md"""
### Implementing Coordinate Descent with normalized features

We seek to obtain a sparse set of weights by minimizing the LASSO cost function

$$\Sigma (prediction - output)^2 + \lambda \times \Sigma_{j=1}^{k} | w_j |$$

(By convention, we do not include $w_0$ in the L1 penalty term. We never want to push the intercept to zero.)

The absolute value sign makes the cost function non-differentiable, so simple gradient descent is not viable (you would need to implement a method called subgradient descent). Instead, we will use **coordinate descent**: at each iteration, we will fix all weights but weight `i` and find the value of weight `i` that minimizes the objective. That is, we look for

$$argmin_{w_i}[\Sigma (prediction - output)^2 + \lambda \times \Sigma_{j=1}^{k} | w_j |]$$

where all weights other than $w_i$ are held to be constant. We will optimize one $w_i$ at a time, circling through the weights multiple times.
  1. Pick a coordinate `i`
  2. Compute $w_i$ that minimizes the cost function $\Sigma[ (prediction - output)^2 ] + \lambda \times \Sigma_{j=1}^{k} | w_j |$
  3. Repeat Steps 1 and 2 for all coordinates, multiple times

For this notebook, we use **cyclical coordinate descent with normalized features**, where we cycle through coordinates 0 to (d-1) in order, and assume the features were normalized as discussed above. The formula for optimizing each coordinate is as follows:

$$w_i = \left\{
    \begin{array}{lr}
      \rho_i + \frac{\lambda}{2} & if \rho_i < -\frac{\lambda}{2} \\
       0  &                        if -\frac{\lambda}{2} \le \rho_i \le \frac{\lambda}{2} \\
      \rho_i - \frac{\lambda}{2} & if \rho_i > \frac{\lambda}{2} \\
    \end{array}
\right.$$

where $\rho_i = \Sigma [ [feature_i] \times (output - prediction + w\_i \times [feature_i]) ]$

Note that we do not regularize the weight of the constant feature (intercept) `w[0]`, so, for this weight, the update is simply:
$w_0 = \rho_i$

"""

# ╔═╡ eb96bfc8-7556-11eb-2045-0f0e7adabe76
md"""
#### Effect of L1 penalty
"""

# ╔═╡ 9c8d9610-754f-11eb-05c3-edd271661b84
begin
	simple_features = [:sqft_living, :bedrooms]
	output_ = :price
	(simple_feature_matrix, output) = get_data(sales, simple_features, output_);
end

# ╔═╡ eb872178-754f-11eb-177d-7bb5ee1046f5
begin
	simple_feature_matrix_n, norms_n = normalize_features(simple_feature_matrix);
end

# ╔═╡ 20799d0c-7550-11eb-16ee-676be88a80af
typeof(simple_feature_matrix), typeof(simple_feature_matrix_n)

# ╔═╡ 06be0f04-7557-11eb-3f7c-bfc5e5a172d4
md"""
We assign some random set of initial weights and inspect the values of ρ[i].

we use `predict_output()` to make predictions on this data.
"""

# ╔═╡ 06a0f408-7557-11eb-182e-8312809145d1
begin
	weights₀ = [1., 4., 1.]

	prediction = predict_output(simple_feature_matrix_n, weights₀);
	prediction, size(prediction)
end

# ╔═╡ 06814074-7557-11eb-2359-8fa03900e3f8
md"""
Compute the values of `ρ[i]` for each feature in this simple model, using the formula given above, using the formula:

$\rho_i = \Sigma [ [feature_i] \times (output - prediction + w\_i \times [feature_i]) ]$

"""

# ╔═╡ fcc38dfc-7557-11eb-048c-93e95ce9bbbc
simple_feature_matrix_n[:, 2]

# ╔═╡ 04b4ddf8-7557-11eb-0d69-f7a89e8936f8
ρ = [
	sum(simple_feature_matrix_n[:, i] .*
		(output - prediction + weights₀[i] * simple_feature_matrix[:, i]))
	for i in 1:length(weights₀)
]

# ╔═╡ 0498db3c-7557-11eb-1af6-9de231a4256f
md"""
**Quiz Question**

Recall that, whenever `ρ[i]` falls between `-l1_penalty/2` and `l1_penalty/2`, the corresponding weight `w[i]` is sent to zero.

Now suppose we were to take one step of coordinate descent on either feature 1 or feature 2.

What range of values of `l1_penalty` **would not** set `w[1]` to zero, but **would** set `w[2]` to zero, if we were to take a step in that coordinate?

"""

# ╔═╡ 0480646e-7557-11eb-3aa9-af64a4502ff9
md"""
we know that by definition:

by definition:
  - for w[2]: $|\rho_2| \le \frac{l_1}{2} \equiv  l_1 \ge 2 * \rho_2$
  - for w[1]: $|\rho_1| \le \frac{l_1}{2} \equiv  l_1 \ge 2 * \rho_1$

"""

# ╔═╡ 046a8ad4-7557-11eb-2aa4-e37c7f37f711
## l1 range that would set w[2] to zero (but not w[1])
min(ρ[2] * 2, ρ[1] * 2)

# ╔═╡ 044fc410-7557-11eb-3a80-eb215ebc3f5f
md"""
**Quiz Question**

What range of values of l1_penalty would set both w[1] and w[2] to zero, if we were to take a step in that coordinate?
"""

# ╔═╡ 0434c2be-7557-11eb-0592-1b733cae76ca
## Any value strictly greater than
max(ρ[2] * 2, ρ[1] * 2)

# ╔═╡ 041c6f52-7557-11eb-050e-c97dee4f3220
md"""
So we can say that `ρ[i]` quantifies the significance of the i-th feature: the larger `ρ[i]` is, the more likely it is for the i-th feature to be retained.
"""

# ╔═╡ 4f9ee984-755a-11eb-06cb-8da12ef68907
md"""
### Single Coordinate Descent Step

Using the formula above, implement coordinate descent that minimizes the cost function over a single feature i. Note that the intercept (weight 0) is not regularized. The function should accept feature matrix, output, current weights, l1 penalty, and index of feature to optimize over. The function should return new weight for feature i.
"""

# ╔═╡ 53782228-755a-11eb-1554-394795823745
function lasso_coordinate_descent_step(ix, f_matrix, output, weights, l₁)
	pred = predict_output(f_matrix, weights)

	ρ_ix = sum(f_matrix[:, ix] .* (output - pred + weights[ix] * f_matrix[:, ix]))

	new_weight_ix = if ix == 1   ## intercept -- do not regularize
		 ρ_ix
	elseif ρ_ix < -l₁ / 2.
		ρ_ix + l₁ / 2.
	elseif ρ_ix >  l₁ / 2.
		ρ_ix - l₁ / 2.
	else
		0.
	end

	new_weight_ix
end

# ╔═╡ 535bf7b0-755a-11eb-06ea-8b0e0baec03c
begin
	ϵ = 1e-7
	exp = -1.7706583789422325
	act = lasso_coordinate_descent_step(
		1,
		[3. / √(13)  1. / √(10); 2. / √(13)  3. / √(10)],
		[1., 1.],
		[1., 4.],
	0.1)

	@test abs(act - exp) ≤ ϵ
end

# ╔═╡ 0403eec8-7557-11eb-18f7-13b0403e0a5b
md"""
#### Cyclical coordinate descent

Now that we have a function that optimizes the cost function over a single coordinate, let us implement cyclical coordinate descent where we optimize coordinates 0, 1, ..., (d-1) in order and repeat.

When do we know to stop? Each time we scan all the coordinates (features) once, we measure the change in weight for each coordinate. If no coordinate changes by more than a specified threshold, we stop.

For each iteration:
1. As you loop over features in order and perform coordinate descent, measure how much each coordinate changes.
2. After the loop, if the maximum change across all coordinates falls below the tolerance, stop. Otherwise, go back to step 1.

Return weights

*IMPORTANT: when computing a new weight for coordinate i, make sure to incorporate the new weights for coordinates 0, 1, ..., i-1. One good way is to update your weights variable in-place.*

"""

# ╔═╡ 03ec653c-7557-11eb-2240-8b060cbb17af
function lasso_cyclical_coordinate_descent(f_matrix, output, init_weights,
		l1_penalty, tolerance)

	weights = copy(init_weights)
	while true
		coord = []

		for ix in 1:length(weights)
			prev_weight_ix = weights[ix]
			weights[ix] = lasso_coordinate_descent_step(ix, f_matrix, output,
				weights, l1_penalty)
			push!(coord, abs(prev_weight_ix - weights[ix]))
		end

		maximum(coord) < tolerance && break
	end
	weights
end

# ╔═╡ dcf8bdc2-755b-11eb-04da-39f89dd5709c
begin
	## utility functions, for the coming questions...

	function features_nzw(features, weights)
  		"""features with non-zero weigths"""
  		return filter(
			((_n, w)=t) -> w > 0.,
			collect(zip([:intercept, features...], weights))
		)

	end

	function features_zw(features, weights, ϵ=1e-7)
  		"""features with zero weigths"""
  		return filter(
			((_n, w)=t) -> w ≤ ϵ,
			collect(zip([:intercept, features...], weights))
			# intercept is constant in the quiz
		)
	end
end

# ╔═╡ 03d065a8-7557-11eb-1cf0-add64e6d7596
begin
	## Using the following parameters, learn the weights on the sales dataset.
	s_features = [:sqft_living, :bedrooms]
	s_outputn = :price
	s_init_weights = zeros(Float64, length(s_features) + 1)  # +1 for intercept
	s_l1_penalty = 1e7
	s_tolerance = 1.0

	s_f_matrix, s_output = get_data(sales, s_features, s_outputn)
	(norm_s_f_matrix, s_norms) = normalize_features(s_f_matrix)

	s_weights = lasso_cyclical_coordinate_descent(norm_s_f_matrix, s_output,
		s_init_weights, s_l1_penalty, s_tolerance)

	s_weights, size(s_weights)
end

# ╔═╡ cf49001e-755e-11eb-33d9-9182f04696f8
md"""
**Quiz Questions**

  1. What is the RSS of the learned model on the normalized dataset? (Hint: use the normalized feature matrix when you make predictions.)
  1. Which features had weight zero at convergence?
"""

# ╔═╡ d35515f6-755e-11eb-3da0-0bc64c7ed3b1
with_terminal() do
	## 1. RSS
	rss = calc_rss(norm_s_f_matrix, s_output, s_weights)
	@printf("rss: %.3f / in scientifc notation: %1.5e\n", rss, rss)
end

# ╔═╡ acd4c536-755e-11eb-1b0c-9573130238dc
## 2. features with weight zero
features_zw(s_features, s_weights)

# ╔═╡ a941a03a-755f-11eb-38f7-455aaa090d7f
md"""
### Evaluating LASSO fit with more features
"""

# ╔═╡ c26ef2fe-755f-11eb-28e4-ef66e958b894
begin
	train_data, test_data = train_test_split(sales; split=0.8, seed=42)
	size(train_data), size(test_data)
end

# ╔═╡ c251d9bc-755f-11eb-2fe2-13c496a0abbd
begin
	a_features = [
		:bedrooms, :bathrooms, :sqft_living, :sqft_lot, :floors,
		:waterfront, :view, :condition, :grade, :sqft_above, :sqft_basement,
		:yr_built, :yr_renovated
	]

	tr_outputn = :price
	(tr_f_matrix, tr_output) = get_data(train_data, a_features, tr_outputn)
	(norm_tr_f_matrix, norms_tr) = normalize_features(tr_f_matrix)

	size(norm_tr_f_matrix), size(norms_tr)
end

# ╔═╡ c0ff7904-755f-11eb-2ad7-572acb587b6f
md"""
##### Learn with l₁ of 1e7
"""

# ╔═╡ c0e5f506-755f-11eb-1ad1-6906e5ccfd85
begin
	l1₁=1e7
	init_weights = zeros(length(a_features) + 1)
	tol = 1.0

	weights1e7 = lasso_cyclical_coordinate_descent(norm_tr_f_matrix, tr_output,
		init_weights, l1₁, tol)
	(size(weights1e7), typeof(weights1e7))
end

# ╔═╡ c0ca4482-755f-11eb-156b-a7011072410a
md"""
**Quiz Question**

What features had non-zero weight in this case?

```julia
features_nzw(a_features, weights1e7)
```
"""

# ╔═╡ 230c74ce-7562-11eb-12d7-43e526f2ec26
with_terminal() do
	for t ∈ features_nzw(a_features, weights1e7)
		@printf("%-20s => %2.5e / %15.5f\n", string(t[1]), t[2], t[2])
	end
end

# ╔═╡ c0b05b74-755f-11eb-319e-9f5eff250eeb
md"""
##### Learn with l₁ of 1e8
"""

# ╔═╡ c095ab50-755f-11eb-2347-79fca41a8555
begin
	l1₂=1e8
	# init_weights2 = zeros(length(a_features) + 1)
	# tol₂ = 1.0

	weights1e8 = lasso_cyclical_coordinate_descent(norm_tr_f_matrix, tr_output,
		init_weights, l1₂, tol)
	(size(weights1e8), typeof(weights1e8))
end

# ╔═╡ c07b70a0-755f-11eb-30ed-cba73d960ff5
md"""
**Quiz Question**

What features had non-zero weight in this case?

```julia
features_nzw(a_features, weights1e8)
```
"""

# ╔═╡ 81964a24-7562-11eb-3ff3-4370b4056f74
with_terminal() do
	for t ∈ features_nzw(a_features, weights1e8)
		@printf("%-20s => %2.5e / %15.5f\n", string(t[1]), t[2], t[2])
	end
end

# ╔═╡ c05c69c6-755f-11eb-0442-ef31019d90ff
md"""
##### Learn with l₁ of 1e4
"""

# ╔═╡ 98a45748-7561-11eb-154a-35b7dd014d41
begin
	l1₃=1e4
	# init_weights₃ = zeros(length(a_features) + 1)
	# tol = 1.0

	weights1e4 = lasso_cyclical_coordinate_descent(norm_tr_f_matrix, tr_output,
		init_weights, l1₃, tol)
	(size(weights1e4), typeof(weights1e4))
end

# ╔═╡ 9856bf2e-7561-11eb-0887-ff4f836e557e
md"""
**Quiz Question**

What features had non-zero weight in this case?

```julia
features_nzw(a_features, weights1e4)
```
"""

# ╔═╡ 08a08158-7563-11eb-0141-15cbef1ee9a6
with_terminal() do
	for t ∈ features_nzw(a_features, weights1e4)
		@printf("%-20s => %2.5e / %15.5f\n", string(t[1]), t[2], t[2])
	end
end

# ╔═╡ a7b5fba4-7563-11eb-2a24-f34b6296b447
md"""
#### Rescaling learned weights

Recall that we normalized our feature matrix, before learning the weights.  To use these weights on a test set, we must normalize the test data in the same way.

Alternatively, we can rescale the learned weights to include the normalization, so we never have to worry about normalizing the test data:

In this case, we must scale the resulting weights so that we can make predictions with *original* features:
 1. Store the norms of the original features to a vector called `norms`:
```
features, norms = normalize_features(features)
```
 2. Run Lasso on the normalized features and obtain a `weights` vector
 3. Compute the weights for the original features by performing element-wise division, i.e.
```
weights_normalized = weights / norms
```
Now, we can apply `weights_normalized` to the test data, without normalizing it!
"""

# ╔═╡ c088e8d0-7563-11eb-0d5c-79ed7f100599
## Create a normalized version of each of the weights learned above. (weights1e4, weights1e7, weights1e8).

begin
all_weights = [
		w ./ vec(norms_tr) for w ∈ (weights1e4, weights1e7, weights1e8)
];

size(all_weights)
end

# ╔═╡ c048f716-7563-11eb-214f-35c773b47c1d
md"""
#### Evaluating each of the learned models on the test data

Let's now evaluate the three models on the test data
"""

# ╔═╡ c02e5654-7563-11eb-390b-ef4352bcf1fa
(test_f_matrix, test_output) = get_data(test_data, a_features, :price);

# ╔═╡ 61567f18-7564-11eb-2301-9dbec94cf157
md"""
Compute the RSS of each of the three normalized weights on the (unnormalized) `test_feature_matrix`
"""

# ╔═╡ 9298e410-7566-11eb-0b86-5376fcf9eabc
md"""
**Quiz Question**

Which model performed best on the test data?
"""

# ╔═╡ 611cc566-7564-11eb-168f-052aa0209973
with_terminal() do
	best_rss, best_weights = (nothing, nothing)
	w_names = (:weights1e4, :weights1e7, :weights1e8)

	for (weights, label) ∈ zip(all_weights, w_names)
  		rss = calc_rss(test_f_matrix, test_output, weights)
  		@printf("for %-10s => rss: %1.5e\n", label, rss)

  		if isnothing(best_rss) || rss < best_rss
    		best_rss, best_weights = rss, label
		end
	end

	@printf("\nbest model/weights: for %-10s => rss: %1.5e\n", best_weights, best_rss)
end

# ╔═╡ 6e03f5ea-7564-11eb-3309-9f6131c552cf

# ╔═╡ Cell order:
# ╟─5b0c59ae-753d-11eb-1a8a-6fa0bf99a357
# ╠═a450c12c-753d-11eb-06fc-29c5d907cb06
# ╟─7e12df34-7549-11eb-30a5-b3099a8f77f6
# ╠═9e5e1060-7549-11eb-0418-55bcd0b8c757
# ╠═a4ce9686-7549-11eb-2172-2d47391d1260
# ╟─2fa16676-754a-11eb-00cc-896d0e28b246
# ╟─4e5d02be-754a-11eb-061f-41ab2f747c54
# ╠═dc99317e-754a-11eb-0a8e-61ec5dff0507
# ╠═e9128af8-7555-11eb-1f4d-d1e17d70f4b3
# ╠═0d7f1fa0-754c-11eb-330b-d31a86f5566c
# ╟─984c778a-7557-11eb-21f8-dd2502fe7b73
# ╟─eb96bfc8-7556-11eb-2045-0f0e7adabe76
# ╠═9c8d9610-754f-11eb-05c3-edd271661b84
# ╠═eb872178-754f-11eb-177d-7bb5ee1046f5
# ╠═20799d0c-7550-11eb-16ee-676be88a80af
# ╟─06be0f04-7557-11eb-3f7c-bfc5e5a172d4
# ╠═06a0f408-7557-11eb-182e-8312809145d1
# ╟─06814074-7557-11eb-2359-8fa03900e3f8
# ╠═fcc38dfc-7557-11eb-048c-93e95ce9bbbc
# ╠═04b4ddf8-7557-11eb-0d69-f7a89e8936f8
# ╟─0498db3c-7557-11eb-1af6-9de231a4256f
# ╟─0480646e-7557-11eb-3aa9-af64a4502ff9
# ╠═046a8ad4-7557-11eb-2aa4-e37c7f37f711
# ╟─044fc410-7557-11eb-3a80-eb215ebc3f5f
# ╠═0434c2be-7557-11eb-0592-1b733cae76ca
# ╟─041c6f52-7557-11eb-050e-c97dee4f3220
# ╟─4f9ee984-755a-11eb-06cb-8da12ef68907
# ╠═53782228-755a-11eb-1554-394795823745
# ╠═535bf7b0-755a-11eb-06ea-8b0e0baec03c
# ╟─0403eec8-7557-11eb-18f7-13b0403e0a5b
# ╠═03ec653c-7557-11eb-2240-8b060cbb17af
# ╠═dcf8bdc2-755b-11eb-04da-39f89dd5709c
# ╠═03d065a8-7557-11eb-1cf0-add64e6d7596
# ╟─cf49001e-755e-11eb-33d9-9182f04696f8
# ╠═d35515f6-755e-11eb-3da0-0bc64c7ed3b1
# ╠═acd4c536-755e-11eb-1b0c-9573130238dc
# ╟─a941a03a-755f-11eb-38f7-455aaa090d7f
# ╠═c26ef2fe-755f-11eb-28e4-ef66e958b894
# ╠═c251d9bc-755f-11eb-2fe2-13c496a0abbd
# ╟─c0ff7904-755f-11eb-2ad7-572acb587b6f
# ╠═c0e5f506-755f-11eb-1ad1-6906e5ccfd85
# ╟─c0ca4482-755f-11eb-156b-a7011072410a
# ╠═230c74ce-7562-11eb-12d7-43e526f2ec26
# ╟─c0b05b74-755f-11eb-319e-9f5eff250eeb
# ╠═c095ab50-755f-11eb-2347-79fca41a8555
# ╟─c07b70a0-755f-11eb-30ed-cba73d960ff5
# ╠═81964a24-7562-11eb-3ff3-4370b4056f74
# ╟─c05c69c6-755f-11eb-0442-ef31019d90ff
# ╠═98a45748-7561-11eb-154a-35b7dd014d41
# ╟─9856bf2e-7561-11eb-0887-ff4f836e557e
# ╠═08a08158-7563-11eb-0141-15cbef1ee9a6
# ╟─a7b5fba4-7563-11eb-2a24-f34b6296b447
# ╠═c088e8d0-7563-11eb-0d5c-79ed7f100599
# ╟─c048f716-7563-11eb-214f-35c773b47c1d
# ╠═c02e5654-7563-11eb-390b-ef4352bcf1fa
# ╟─61567f18-7564-11eb-2301-9dbec94cf157
# ╟─9298e410-7566-11eb-0b86-5376fcf9eabc
# ╠═611cc566-7564-11eb-168f-052aa0209973
# ╠═6e03f5ea-7564-11eb-3309-9f6131c552cf


text:
using BinDeps
using Compat

# Libdl has been moved to a standard library module
using Compat.Libdl

@BinDeps.setup

# Add library dependency for direct method
osqp = library_dependency("osqp", aliases=["libosqp"])

# Current version
version = "0.4.0"
qdldl_version = "0.1.2"

# Get current operating system
osqp_platform =
if Compat.Sys.islinux()
    "linux"
elseif Compat.Sys.isapple()
    "mac"
elseif Compat.Sys.iswindows()
    "windows"
else
    error("Platform not supported!")
end

# Provide binaries for each operating system
archive_name="osqp-$version-$osqp_platform$(Sys.WORD_SIZE)"
provides(Binaries, URI("https://dl.bintray.com/bstellato/generic/OSQP/$version/$archive_name.tar.gz"), [osqp], unpacked_dir="$archive_name/lib", os=:Darwin)
provides(Binaries, URI("https://dl.bintray.com/bstellato/generic/OSQP/$version/$archive_name.tar.gz"), [osqp], unpacked_dir="$archive_name/lib", os=:Windows)

# Build from sources on Linux
provides(Sources, URI("https://github.com/oxfordcontrol/osqp/archive/v$version.tar.gz"),
    [osqp], unpacked_dir="osqp-$version", os = :Linux)

# Define directories locations deps/usr and deps/src
prefix = joinpath(BinDeps.depsdir(osqp),"usr")
srcdir = joinpath(BinDeps.depsdir(osqp),"src","osqp-$version")
blddir = joinpath(srcdir, "build")
downloadsdir = joinpath(BinDeps.depsdir(osqp),"downloads")

# Define library name
libname = "libosqp.$(dlext)"

# File move rule to move QDLDL sources
struct FileMoveRule <: BinDeps.BuildStep
    src::AbstractString
    dest::AbstractString
end
Base.run(fc::FileMoveRule) = isfile(fc.dest) || mv(fc.src, fc.dest, remove_destination=true)

provides(SimpleBuild,
    (@build_steps begin
    GetSources(osqp)
    # Copy QDLDL sources inside OSQP sources
    FileDownloader("https://github.com/oxfordcontrol/qdldl/archive/v$qdldl_version.tar.gz",
		   joinpath(downloadsdir, "external", "qdldl-$qdldl_version.tar.gz"))
    FileUnpacker(joinpath(downloadsdir, "external", "qdldl-$qdldl_version.tar.gz"),
		 joinpath(srcdir, downloadsdir),
		 "qdldl-$qdldl_version")
    FileMoveRule(joinpath(srcdir, downloadsdir, "qdldl-$qdldl_version"),
		 joinpath(srcdir, "lin_sys", "direct", "qdldl", "qdldl_sources"))
    CreateDirectory(joinpath(prefix, "lib"))
    @build_steps begin
        CreateDirectory(blddir)
        @build_steps begin
        ChangeDirectory(blddir)
        FileRule(joinpath(prefix, "lib", libname), @build_steps begin
            `cmake -G "Unix Makefiles" ..`
            `make osqp`
            `mv out/$libname $prefix/lib`
        end)
        end
    end
    end),
[osqp], os = :Linux)

@BinDeps.install Dict(:osqp => :osqp)


text:
"""
This module provides ways to parse and write files in
[FASTA format](http://en.wikipedia.org/wiki/FASTA_format) in Julia.
It is designed to be lightweight and fast; the parsing method is inspired by
[kseq.h](http://lh3lh3.users.sourceforge.net/kseq.shtml). It can read and write
files on the fly, keeping only one entry at a time in memory, and it can read and
write gzip-compressed files.

See [`FastaReader`](@ref), [`FastaWriter`](@ref), [`readfasta`](@ref),
[`writefasta`](@ref), [`readentry`](@ref) and [`writeentry`](@ref).
"""
module FastaIO

using GZip

export
    FastaReader,
    readentry,
    rewind,
    readfasta,
    FastaWriter,
    writeentry,
    writefasta

import Base.close, Base.show, Base.eof, Base.write
import Base: iterate

const fasta_buffer_size = 4096

mutable struct FastaReader{T}
    # public but read-only
    num_parsed::Int        # number of parsed entries so far
    # private
    f::IO
    is_eof::Bool           # did we reach end of file?
    rbuffer::Vector{UInt8} # read buffer
    rbuf_sz::Int           # read buffer size
    rbuf_pos::Int          # read buffer cursor
    lbuffer::Vector{UInt8} # line buffer
    lbuf_sz::Int           # line buffer size
    mbuffer::Vector{UInt8} # multi-line buffer
    mbuf_sz::Int           # multi-line buffer size
    own_f::Bool
    function FastaReader{T}(filename::AbstractString) where {T}
        fr = new{T}(0, gzopen(filename), false,
                    Array{UInt8}(undef, fasta_buffer_size), 0, 0,
                    Array{UInt8}(undef, fasta_buffer_size), 0,
                    Array{UInt8}(undef, fasta_buffer_size), 0,
                    true)
        finalizer(close, fr)
        return fr
    end
    function FastaReader{T}(io::IO) where {T}
        new{T}(0, io, false,
               Array{UInt8}(undef, fasta_buffer_size), 0, 0,
               Array{UInt8}(undef, fasta_buffer_size), 0,
               Array{UInt8}(undef, fasta_buffer_size), 0,
               false)
    end
end

"""
    FastaReader{T}(file::Union{AbstractString,IO})

This creates an object which is able to parse FASTA files, one entry at a time. `file` can be a plain text
file or a gzip-compressed file (it will be autodetected from the content). The type `T` determines the output
type of the sequences (see [The sequence storage type](@ref) section for more information) and it defaults to
`String`.

The data can be read out by iterating the `FastaReader` object:

```julia
for (name, seq) in FastaReader("somefile.fasta")
    # do something with name and seq
end
```

As shown, the iterator returns a tuple containing the description (always a `String`) and the data (whose
type is set when creating the `FastaReader` object (e.g. `FastaReader{Vector{UInt8}}(filename)`).

The `FastaReader` type has a field `num_parsed` which contains the number of entries parsed so far.

Other ways to read out the data are via the [`readentry`](@ref) and [`readfasta`](@ref) functions.
"""
FastaReader(file::Union{AbstractString,IO}) = FastaReader{String}(file)

"""
    FastaReader(f::Function, filename::AbstractString, [sequence_type::Type = String])

This format of the constructor is useful for do-notation, i.e.:

```julia
FastaReader(filename) do fr
    # read out the data from fr, e.g.
    for (name, seq) in fr
        # do something with name and seq
    end
end
```

which ensures that the [`close`](@ref) function is called and is thus recommended (otherwise
the file is closed by the garbage collector when the `FastaReader` object goes out of scope).
"""
function FastaReader(f::Function, filename::AbstractString, T::Type=String)
    fr = FastaReader{T}(filename)
    try
        f(fr)
    finally
        close(fr)
    end
end

"""
    close(fr::FastaReader)

This function extends `Base.close` and closes the stream associated with the `FastaReader`; the reader must not
be used any more after this function is called.
"""
close(fr::FastaReader) = fr.own_f && close(fr.f)

"""
    rewind(fr::FastaReader)

This function rewinds the reader, so that it can restart the parsing again without closing and re-opening it.
It also resets the value of the `num_parsed` field.
"""
function rewind(fr::FastaReader)
    seek(fr.f, 0)
    fr.is_eof = false
    fr.num_parsed = 0
    fr.rbuf_sz = 0
    fr.rbuf_pos = 0
    fr.lbuf_sz = 0
    fr.mbuf_sz = 0
    return
end

read_chunk_ll(fr::FastaReader, s::GZipStream) = gzread(s, pointer(fr.rbuffer), fasta_buffer_size)
read_chunk_ll(fr::FastaReader, s::IOStream) =
    ccall(:ios_readall, UInt, (Ptr{Nothing}, Ptr{Nothing}, UInt), fr.f.ios, fr.rbuffer, fasta_buffer_size)
function read_chunk_ll(fr::FastaReader, s::IO)
    ret = 0
    while !eof(fr.f) && ret < fasta_buffer_size
        ret += 1
        fr.rbuffer[ret] = read(fr.f, UInt8)
    end
    return ret
end

function read_chunk(fr::FastaReader)
    if fr.is_eof
        return
    end
    ret = read_chunk_ll(fr, fr.f)
    ret == -1 && error("read failure")
    fr.rbuf_sz = ret
    fr.rbuf_pos = 1
    if ret == 0
        fr.is_eof = true
    end
    return
end

function readline(fr::FastaReader)
    fr.lbuf_sz = 0
    found = false
    while !fr.is_eof
        if fr.rbuf_pos == 0
            read_chunk(fr::FastaReader)
        end
        i = fr.rbuf_pos
        cr = false
        while i <= fr.rbuf_sz
            c = fr.rbuffer[i]
            if c == UInt8('\n')
                found = true
                break
            else
                cr = (c == UInt8('\r'))
            end
            i += 1
        end
        i -= 1 + cr
        chunk_len = i - fr.rbuf_pos + 1
        free_sbuf = length(fr.lbuffer) - fr.lbuf_sz
        gap = chunk_len - free_sbuf
        if gap > 0
            resize!(fr.lbuffer, length(fr.lbuffer) + gap)
        end

        #fr.lbuffer[fr.lbuf_sz + (1:chunk_len)] = fr.rbuffer[fr.rbuf_pos:i]
        copyto!(fr.lbuffer, fr.lbuf_sz + 1, fr.rbuffer, fr.rbuf_pos, chunk_len)
        fr.lbuf_sz += chunk_len

        i += 2 + cr
        if i > fr.rbuf_sz
            i = 0
        end
        fr.rbuf_pos = i
        found && break
    end
    if !fr.is_eof && fr.lbuf_sz == 0
        # empty line, ignore and run again
        return readline(fr)
    end
    return
end

# gets the ID from lbuffer and reads FASTA entry into mbuffer
function _next_step(fr::FastaReader)
    if fr.lbuffer[1] != UInt8('>')
        error("invalid FASTA file: description does not start with '>'")
    end
    if fr.lbuf_sz == 1
        error("invalid FASTA file: empty description")
    end
    name = String(fr.lbuffer[2:fr.lbuf_sz])
    isascii(name) || error("invalid non-ASCII description in FASTA file")
    fr.mbuf_sz = 0
    while true
        readline(fr)
        if fr.lbuf_sz == 0 || fr.lbuffer[1] == UInt8('>')
            break
        end
        gap = fr.lbuf_sz - (length(fr.mbuffer) - fr.mbuf_sz)
        if gap > 0
            resize!(fr.mbuffer, length(fr.mbuffer) + gap)
        end
        #fr.mbuffer[fr.mbuf_sz + (1:fr.lbuf_sz)] = fr.lbuffer[1:fr.lbuf_sz]
        copyto!(fr.mbuffer, fr.mbuf_sz + 1, fr.lbuffer, 1, fr.lbuf_sz)
        fr.mbuf_sz += fr.lbuf_sz
    end
    return name
end

# extracts the sequence of the FASTA entry from the mbuffer, the result is of type T
_next_seq(fr::FastaReader{Vector{UInt8}}) = fr.mbuffer[1:fr.mbuf_sz]
_next_seq(fr::FastaReader{String}) = ccall(:jl_pchar_to_string, Ref{String}, (Ptr{UInt8},Int), fr.mbuffer, fr.mbuf_sz)
_next_seq(fr::FastaReader{T}) where T = T(fr.mbuffer[1:fr.mbuf_sz])

function _next(fr::FastaReader{T}) where T
    name = _next_step(fr)
    seq = _next_seq(fr)::T
    fr.num_parsed += 1
    return (name, seq)
end

Base.eltype(fr::FastaReader{T}) where T = Tuple{String, T}
Base.IteratorSize(fr::FastaReader) = Base.SizeUnknown()

function iterate(fr::FastaReader)
    rewind(fr)
    readline(fr)
    if fr.lbuf_sz == 0
        error("empty FASTA file")
    end
    return fr.is_eof ? nothing : (_next(fr), nothing)
end
function iterate(fr::FastaReader, x::Nothing)
    fr.is_eof && return nothing
    return _next(fr), nothing
end

"""
    readentry(fr::FastaReader)

This function can be used to read entries one at a time:

```julia
fr = FastaReader("somefile.fasta")
name, seq = readentry(fr)
```

See also the [`eof`](@ref) function.
"""
function readentry(fr::FastaReader)
    fr.is_eof && throw(EOFError())
    if fr.num_parsed == 0
        readline(fr)
        if fr.lbuf_sz == 0
            error("empty FASTA file")
        end
    end
    item = _next(fr)
    return item
end

"""
    eof(fr::FastaReader)

This function extends `Base.eof` and tests for end-of-file condition; it is useful when using
[`readentry`](@ref):

```julia
fr = FastaReader("somefile.fasta")
while !eof(fr)
   name, seq = readentry(fr)
   # do something
end
close(fr)
```
"""
eof(fr::FastaReader) = fr.is_eof

function show(io::IO, fr::FastaReader{T}) where {T}
    print(io, "FastaReader(input=\"$(fr.f)\", out_type=$T, num_parsed=$(fr.num_parsed), eof=$(fr.is_eof))")
end

"""
    readfasta(file::Union{String,IO}, [sequence_type::Type = String])

This function parses a whole FASTA file at once and stores it into memory. The result is a `Vector{Any}`
whose elements are tuples consisting of `(description, sequence)`, where `description` is a
`String` and `sequence` contains the sequence data, stored in a container type defined by
the `sequence_type` optional argument (see [The sequence storage type](@ref) section for more information).
"""
readfasta(filename::AbstractString, ::Type{T}=String) where T =
    gzopen(io -> readfasta(io, T), filename)
readfasta(io::IO, ::Type{T}=String) where T = collect(FastaReader{T}(io))

mutable struct FastaWriter
    f::IO
    in_seq::Bool
    entry_chars::Int
    desc_chars::Int
    parsed_nl::Bool
    pos::Int
    entry::Int
    own_f::Bool
    at_start::Bool
    function FastaWriter(io::IO)
        fw = new(io, false, 0, 0, false, 0, 1, false, true)
        finalizer(close, fw)
        return fw
    end
    function FastaWriter(filename::AbstractString, mode::AbstractString = "w")
        if endswith(filename, ".gz")
            of = gzopen
        else
            of = open
        end
        fw = new(of(filename, mode), false, 0, 0, false, 0, 1, true, true)
        finalizer(close, fw)
        return fw
    end
end

FastaWriter() = FastaWriter(STDOUT)

"""
    FastaWriter(filename::AbstractString, [mode::String = "w"])
    FastaWriter([io::IO = STDOUT])
    FastaWriter(f::Function, args...)

This creates an object which is able to write formatted FASTA files which conform to the specifications
detailed in the section titled [The FASTA format](@ref), via the [`write`](@ref) and [`writeentry`](@ref)
functions.

The third form allows to use do-notation:

```julia
FastaWriter("somefile.fasta") do fw
    # write the file
end
```

which is strongly recommended since it ensures that the [`close`](@ref) function is called at the end of
writing: this is crucial, as failing to do so may result in incomplete files (this is done by the
finalizer, so it will still happen automatically if the `FastaWriter` object goes out of scope and is
garbage-collected, but there is no guarantee that this will happen if Julia exits).

If the `filename` ends with `.gz`, the result will be gzip-compressed.

The `mode` flag can be used to set the opening mode of the file; use `"a"` to append to an existing
file.

The `FastaWriter` object has an `entry::Int` field which stores the number of the entry which is
currently being written.
"""
function FastaWriter(f::Function, args...)
    fw = FastaWriter(args...)
    try
        f(fw)
    finally
        close(fw)
    end
end

"""
    write(fw::FastaWriter, item)

This function extends `Base.write` and streams items to a FASTA file, which will be formatted
according to the specifications detailed in the section titled [The FASTA format](@ref).

When using this method, description lines are marked by the fact that they begin with a `'>'`
character; anything else is assumed to be part of the sequence data.

If `item` is a `Vector`, `write` will be called iteratively over it; if it is a `String`,
a newline will be appended to it and it will be dumped. For example the following code:

```julia
FastaWriter("somefile.fasta") do fw
    for s in [">GENE1", "GCA", "TTT", ">GENE2", "ATTAGC"]
        write(fw, s)
    end
end
```

will result in the file:

```text
>GENE1
GCATTT
>GENE2
ATTAGC
```

If `item` is not a `Vector` nor a `String`, it must be convertible to an ASCII character, and
it will be piped into the file. For example the following code:

```julia
data = \"\"\"
  >GENE1
  GCA
  TTT
  >GENE2
  ATT
  AGC
  \"\"\"

FastaWriter("somefile.fasta") do fw
    for ch in data
        write(fw, ch)
    end
end
```

will result in the same file as above.
"""
function write(fw::FastaWriter, c)
    ch = convert(Char, c)
    isascii(ch) || error("invalid (non-ASCII) character: $c (entry $(fw.entry) of FASTA input)")
    if ch == '\n' && !fw.at_start
        fw.parsed_nl = true
        if !fw.in_seq
            fw.desc_chars == 1 && error("empty description (entry $(fw.entry) of FASTA input")
            write(fw.f, '\n')
            fw.pos = 0
            fw.in_seq = true
        end
    end
    isspace(ch) && (fw.at_start || fw.in_seq || fw.desc_chars <= 1) && return
    fw.at_start && ch != '>' && error("no desctiption given (entry $(fw.entry) of FASTA input")
    fw.at_start = false
    if fw.parsed_nl
        @assert fw.in_seq
        if ch == '>'
            fw.entry_chars > 0 || error("description must span a single line (entry $(fw.entry) of FASTA input)")
            write(fw.f, '\n')
            fw.in_seq = false
            fw.pos = 0
            fw.entry += 1
            fw.entry_chars = 0
            fw.desc_chars = 0
        end
    elseif fw.in_seq && ch == '>'
        error("character '>' not allowed in sequence data (entry $(fw.entry) of FASTA input)")
    end
    if fw.pos == 80
        if !fw.in_seq
            warn("description line longer than 80 characters (entry $(fw.entry) of FASTA input)")
        else
            write(fw.f, '\n')
            fw.pos = 0
        end
    end
    write(fw.f, ch)
    fw.pos += 1
    if fw.in_seq
        fw.entry_chars += 1
    else
        fw.desc_chars += 1
    end
    fw.parsed_nl = false
    return
end

function write(fw::FastaWriter, s::Vector)
    for c in s
        write(fw, c)
    end
end

function write(fw::FastaWriter, s::AbstractString)
    for c in s
        write(fw, c)
    end
    write(fw, '\n')
end

"""
    writeentry(fw::FastaWriter, description::AbstractString, sequence)

This function writes one entry to the FASTA file, following the specifications detailed in
the section titled [The FASTA format](@ref). The `description` is without the initial `'>'` character.
The `sequence` can be any iterable object whose elements are convertible to ASCII characters.

Example:

```julia
FastaWriter("somefile.fasta") do fw
    for (desc,seq) in [("GENE1", "GCATT"), ("GENE2", "ATTAGC")]
        writeentry(fw, desc, seq)
    end
end
```
"""
function writeentry(fw::FastaWriter, desc::AbstractString, seq)
    !fw.at_start && write(fw, '\n')
    desc = strip(String(desc))
    isascii(desc) || error("description must be ASCII (entry $(fw.entry+1) of FASTA input)")
    if findfirst(isequal('\n'), desc) ∉ (nothing, 0) # TODO: remove 0 when support for julia 0.6 is dropped
        error("newlines are not allowed within description (entry $(fw.entry+1) of FASTA input)")
    end
    write(fw, '>')
    write(fw, desc)
    write(fw, '\n')
    #write(fw, seq)
    #write(fw, '\n')
    fw.entry_chars = writefastaseq(fw.f, seq, fw.entry, false)
    fw.in_seq = true
    fw.parsed_nl = false
    fw.pos = 0
    fw.entry_chars > 0 || error("empty sequence data (entry $(fw.entry) of FASTA input)")
    return
end

"""
    close(fw::FastaWriter)

This function extends `Base.close` and it should always be explicitly used for finalizing
the `FastaWriter` once the writing has finished, unless the do-notation is used when creating
it.
"""
function close(fw::FastaWriter)
    try
        write(fw.f, '\n')
        flush(fw.f)
    catch err
        isa(err, EOFError) || rethrow(err)
    end
    fw.pos = 0
    fw.parsed_nl = true
    fw.own_f && close(fw.f)
    return
end

function show(io::IO, fw::FastaWriter)
    print(io, "FastaWriter(input=\"$(fw.f)\", entry=$(fw.entry)")
end

function writefastaseq(io::IO, seq, entry::Int, nl::Bool = true)
    i = 0
    entry_chars = 0
    for c in seq
        if i == 80
            write(io, '\n')
            i = 0
        end
        ch = convert(Char, c)
        isascii(ch) || error("invalid (non-ASCII) character: $c (entry $entry of FASTA input)")
        isspace(ch) && continue
        ch != '>' || error("character '>' not allowed in sequence data (entry $entry of FASTA input)")
        write(io, ch)
        i += 1
        entry_chars += 1
    end
    nl && write(io, '\n')
    return entry_chars
end

"""
    writefasta([io::IO = STDOUT], data)

This version of the function writes to an already opened `IO` stream, defaulting to `STDOUT`.
"""
function writefasta(io::IO, data)
    entry = 0
    for (desc, seq) in data
        entry += 1
        desc = strip(String(desc))
        isascii(desc) || error("description must be ASCII (entry $entry of FASTA input)")
        if isempty(desc)
            error("empty description (entry $entry of FASTA input")
        end
        if findfirst(isequal('\n'), desc) ∉ (nothing, 0) # TODO: remove 0 when support for julia 0.6 is dropped
            error("newlines are not allowed within description (entry $entry of FASTA input)")
        end
        if length(desc) > 79
            warn("description line longer than 80 characters (entry $entry of FASTA input)")
        end
        println(io, ">", desc)
        entry_chars = writefastaseq(io, seq, entry)
        entry_chars > 0 || error("empty sequence data (entry $entry of FASTA input)")
    end
end
writefasta(data) = writefasta(STDOUT, data)

"""
    writefasta(filename::String, data, [mode::String = "w"])

This function dumps data to a FASTA file, auto-formatting it so to follow the specifications detailed in
the section titled [The FASTA format](@ref). The `data` can be anything which is iterable and which produces
`(description, sequence)` tuples upon iteration, where the `description` must be convertible to
a `String` and the `sequence` can be any iterable object which yields elements convertible
to ASCII characters (e.g. a `String`, a `Vector{UInt8}` etc.).

Examples:

```julia
writefasta("somefile.fasta", [("GENE1", "GCATT"), ("GENE2", "ATTAGC")])
writefasta("somefile.fasta", ["GENE1" => "GCATT", "GENE2" => "ATTAGC"])
```

If the `filename` ends with `.gz`, the result will be a gzip-compressed file.

The `mode` flag determines how the `filename` is open; use `"a"` to append the data to an existing
file.
"""
function writefasta(filename::AbstractString, data, mode::AbstractString = "w")
    if endswith(filename, ".gz")
        of = gzopen
    else
        of = open
    end
    of(filename, mode) do f
        writefasta(f, data)
    end
end

end


text:
const nltk_atoms = collect.(["--", "...", "``", "\$"])
const nltk_suffixes = collect.(["'ll", "'re", "'ve", "n't", "'s", "'m", "'d"])
const nltk_splits = [
    "cannot"=>("can", "not"),
    "gimme"=>("gim", "me"),
    "lemme"=>("lem", "me"),
    "mor'n"=>("mor", "'n"),
    "d'ye"=>("d'n", "e"),
    "gonna"=>("gon","na"),
    "gotta"=>("got","ta"),
    "wanna"=>("wan","na"),
    "'tis"=>("'t", "is"),
    "'twas"=>("'t", "was"),
]

"""
    nltk_word_tokenize(input)

NLTK's word tokenizer.
It is an extension on the Punctuation Preserving Penn Treebank tokenizer,
mostly to better handle unicode.

Punctuation is still preserved as its own token.
This includes periods which will be stripped from words.

The tokenizer still seperates out contractions:
"shouldn't" becomes ["should", "n't"]

The input should be a single sentence;
but again it will likely be relatively fine if it isn't.
Depends exactly what you want it for.

This matches to the most commonly used `nltk.word_tokenize`, minus the sentence tokenizing step.
"""
function nltk_word_tokenize(input)
    ts = TokenBuffer(input)
    isempty(input) && return ts.tokens
    stop = ts.input[end] == '.' # `.` is usually absorbed into tokens (`Dr.`)
                              # Treat the last `.` specially.
    stop && pop!(ts.input)
    while !isdone(ts)
        spaces(ts) && continue
        openquote(ts) ||
        suffixes(ts, nltk_suffixes) ||
        atoms(ts, nltk_atoms) ||
        replaces(ts, nltk_splits) ||
        number(ts) ||
        character(ts)
        !isdone(ts) && closingquote(ts)
    end
    stop && push!(ts.tokens, ".")
    return ts.tokens
end


text:
include("imgtransform.jl")
using .ImgTransform
using ImageTransformations
using Images
using DelimitedFiles
using Dates
using Colors

randStr(len::Int64) = String(UInt8.(rand([97:122; 65:90; 48:57], len)))

log(msg::String) = println("$(Dates.now()) $msg")

function processInput(inputDir::String, transformations, outputDir::String)
    paths = joinpath.(inputDir, readdir(inputDir))

    for path in paths
        log("Loading image: $path")
        image = load(path)
        log("Binarization")
        image = ImgTransform.imgToBitMatrix(image)
        log("Cutting out shape")
        image = ImgTransform.cutoutShape(image)
        image = xor.(image, 1)

        for (i, transform) in enumerate(transformations)
            log("Selecting $i-th transformation")
           # log("Looking for available random name")
            newpath = joinpath(outputDir, randStr(rand(5:10))*".bmp")
            #while isfile(path)
            #    newpath = joinpath(outputDir, randStr(rand(1:15)))
            #end

            log("Transforming")
            transformed = transform(image)
            result = Gray.(transformed)

            log("Writing transformed image to: $newpath")
            save(newpath, result)
        end
    end
end

inputs = ["/home/jakub/Dokumenty/skola/siete/neural/presets/circles/", "/home/jakub/Dokumenty/skola/siete/neural/presets/rectangles/", "/home/jakub/Dokumenty/skola/siete/neural/presets/triangles/"] # input directories
outputs = ["/run/media/jakub/3315-3BCD/trainset3/circles", "/run/media/jakub/3315-3BCD/trainset3/rectangles", "/run/media/jakub/3315-3BCD/trainset3/triangles"] # output directories

rotations = [x -> imrotate(x, θ, fillvalue=1) for θ in (10*π/180):(10*π/180):2π]
resize = [i -> imresize(i, 30, 30)]
rotresize = map(rot -> x -> imresize(rot(x), 30, 30), rotations)
#TODO: add more transfomations; gaussian noise, add blob to random place, remove random part of the shape

transformations = [resize]
#transformations = [rotations]
#transformations = [[x -> imrotate(x, π/4, fillvalue=1)]]

log("Creating output directories")
mkpath.(outputs)

for (input, output) in zip(inputs, outputs)
    log("Entering directory: $input")
    for (i, transforms) in enumerate(transformations)
        log("\tUsing transformations set: $i")
        processInput(input, transforms, output)
    end
end

max_stars_repo_name: achjaj/shape-recognition

max_issues_repo_name: achjaj/shape-recognition

max_forks_repo_name: achjaj/shape-recognition


text:
using Documenter
using TexTables

makedocs(
    modules = [TexTables],
    sitename= "TexTables.jl",
    format  = :html,
    clean   = false,
    authors = "Jacob Adenbaum",
    pages   = [
              "Introduction"=>  "index.md",
              "Easy Examples"=> "easy_examples.md",
              "Basic Usage" =>  "basic_usage.md",
              "Regression API"=> "regression_tables.md",
              "Advanced Usage"=>"advanced_usage.md"]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo   = "github.com/jacobadenbaum/TexTables.jl.git",
    julia  = "0.6",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
function vtk_grid(dtype::VTKImageData, filename::AbstractString,
                  Nx::Integer, Ny::Integer, Nz::Integer=1;
                  origin = (0.0, 0.0, 0.0),
                  spacing = (1.0, 1.0, 1.0),
                  extent = nothing, kwargs...)
    Npts = Nx*Ny*Nz
    Ncls = num_cells_structured(Nx, Ny, Nz)
    ext = extent_attribute(Nx, Ny, Nz, extent)

    xvtk = XMLDocument()
    vtk = DatasetFile(dtype, xvtk, filename, Npts, Ncls; kwargs...)

    # VTKFile node
    xroot = vtk_xml_write_header(vtk)

    # ImageData node
    xGrid = new_child(xroot, vtk.grid_type)
    set_attribute(xGrid, "WholeExtent", ext)

    No = length(origin)
    Ns = length(spacing)

    if No > 3
        throw(ArgumentError("origin array must have length <= 3"))
    elseif Ns > 3
        throw(ArgumentError("spacing array must have length <= 3"))
    end

    origin_str = join(origin, " ")
    spacing_str = join(spacing, " ")

    while No != 3
        # Fill additional dimensions (e.g. the z dimension if 2D grid)
        origin_str *= (" 0.0")
        No += 1
    end

    while Ns != 3
        spacing_str *= (" 1.0")
        Ns += 1
    end

    set_attribute(xGrid, "Origin", origin_str)
    set_attribute(xGrid, "Spacing", spacing_str)

    # Piece node
    xPiece = new_child(xGrid, "Piece")
    set_attribute(xPiece, "Extent", ext)

    vtk
end

vtk_grid(filename::AbstractString, xyz::Vararg{Integer}; kwargs...) =
    vtk_grid(VTKImageData(), filename, xyz...; kwargs...)

"""
    vtk_grid(filename, x::AbstractRange{T}, y::AbstractRange{T}, [z::AbstractRange{T}];
             kwargs...)

Create image data (`.vti`) file.

Along each direction, the grid is specified in terms of an AbstractRange object.

# Examples

```jldoctest
julia> vtk = vtk_grid("abc", 1:0.2:5, 2:1.:3, 4:1.:5)  # 3D dataset
VTK file 'abc.vti' (ImageData file, open)

julia> vtk = vtk_grid("abc", 1:5, 2:1.:3, range(4, 5; length = 3))  # different kinds of ranges
VTK file 'abc.vti' (ImageData file, open)

julia> vtk = vtk_grid("abc", 1:0.2:5, 2:1.:3)  # 2D dataset
VTK file 'abc.vti' (ImageData file, open)

julia> vtk = vtk_grid("def",
                      LinRange(0., 5., 10),
                      LinRange(0., 2π, 16),
                      LinRange(1., 10., 12))
VTK file 'def.vti' (ImageData file, open)

```

"""
function vtk_grid(filename::AbstractString, xyz::Vararg{AbstractRange}; kwargs...)
    Nxyz = length.(xyz)
    origin = first.(xyz)
    spacing = step.(xyz)
    vtk_grid(VTKImageData(), filename, Nxyz...;
             origin=origin, spacing=spacing, kwargs...)
end


text:
using DifferentialEquations
using LinearAlgebra
using Plots
using DeltaRobot

# u0 =  [pi/4; pi/2; 0.0; pi/4; pi/2; 0.0; pi/4; pi/2; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]; # what should lambda be?
u0 =  [pi/4; 3*pi/4; 0.0; pi/4; 3*pi/4; 0.0; pi/4; 3*pi/4; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]; # what should lambda be?
du0 = [0.0; 0.0 ; 0.0 ; 0.0; 0.0 ; 0.0 ; 0.0; 0.0 ; 0.0 ; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]; # why must I give this?
differential_vars = [true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, false, false, false, false, false, false]
@assert length(u0) == length(du0)
@assert length(u0) == length(differential_vars)
@assert u0[10:18] == du0[1:9] # make sure initial conditions are consistent
p = DeltaRobotArmParams(0.3, 0.2, 0.1, 0.5, 0.5, 9.8, 0.09, 0.09, 0.09);
tspan = (0.0, 20.0)
prob = DAEProblem(delta_robot!, du0, u0, tspan, p, differential_vars=differential_vars)
dt = 0.05
sol = solve(prob, saveat=dt)

gr()
anim = animate_robot(sol, p)
gif(anim, "animations/robot_1.gif", fps=convert(Int64, 1 / dt))

plotly()
plot(sol, vars=(1, 2, 3))
plot(sol.t, sol[1,:])
plot(sol.t, sol[2,:])
plot(sol.t, sol[3,:])

gr()
prob_damped = DAEProblem(delta_robot_damped!, du0, u0, tspan, p, differential_vars=differential_vars)
sol_damped = solve(prob_damped, saveat=dt)
anim = animate_robot(sol_damped, p)
gif(anim, "animations/damp_robot_1.gif", fps=convert(Int64, 1 / dt))

max_stars_repo_name: zborffs/Delta

max_issues_repo_name: zborffs/Delta

max_forks_repo_name: zborffs/Delta


text:
ENV["DYLD_LIBRARY_PATH"]="/usr/lib/"

using ConvDiffMIPDECO
# using PyPlot
using jInv.Mesh
using jInv.ForwardShare
using jInv.InverseSolve
# using jInvVis
using jInv.LinearSolvers
using MAT
using MUMPSjInv

# filename= "2DmodelLShaped.mat"
filename= joinpath(dirname(pathof(ConvDiffMIPDECO)),"..","examples","Peaks2D.mat")
file = matread(filename)
wtrue = vec(file["W"]) ## Sources on finest mesh
domain = file["domain"]
m      = file["m"]

M      = getRegularMesh(domain,m)
# viewImage2D(wtrue,M)

rec = rand(200,2) .*[domain[2] domain[4]]
# plot(rec[:,1],rec[:,2],".r")

# build linear interpolation matrix from nodes to receiver locations
x1,x2 = getNodalAxes(M)
P = interpmat(x1,x2, rec);

# build param that holds forward problem
pFor = getConvDiffFEMParam(M,Ainv=getMUMPSsolver(),P=P);

dtrue, pFor = getData(vec(wtrue),pFor);
utrue = pFor.Fields

file["utrue"]  = reshape(utrue,M.n[1]+1,M.n[2]+1);
file["dtrue"]  = dtrue
file["rec"]    = rec
file["sig"]    = pFor.sig
file["v"]      = pFor.v

matwrite(filename,file)

# viewImage2D(utrue,getPaddedMesh(M))
#plot(rec[:,1],rec[:,2],".r")
#
# title("true fields")

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
const testdatadir = joinpath(dirname(@__FILE__), "data")

REPO_URL = "https://github.com/yeesian/ArchGDALDatasets/blob/master/"

remotefiles = [
    "data/A.tif",
    "data/utmsmall.tif",
    "gdalworkshop/world.tif",
    "ospy/data4/aster.img",
    # "ospy/data4/aster.rrd",
    "ospy/data5/aster.img",
    # "ospy/data5/aster.rrd",
    "ospy/data5/doq1.img",
    # "ospy/data5/doq1.rrd",
    "ospy/data5/doq2.img",
    # "ospy/data5/doq2.rrd",
    "pyrasterio/example.tif",
    "pyrasterio/example2.tif",
    "pyrasterio/example3.tif",
    "pyrasterio/float_nan.tif",
    "pyrasterio/float.tif",
    "pyrasterio/RGB.byte.tif",
    "pyrasterio/shade.tif"
]

for f in remotefiles
    # create the directories if they don't exist
    currfile = joinpath(testdatadir, f)
    currdir = dirname(currfile)
    isdir(currdir) || mkpath(currdir)
    # download the file
    isfile(currfile) || download(REPO_URL * f * "?raw=true", currfile)
end


text:
    """
there always has been a source of fascination
with a natural Network model - abstract here
used to brainstorm new Ideas
    """

    struct NaturalNet <: NN #extends the most abstract model in NeuralNet directory
    synapses::Any #  are always between the axon of one neuron and the soma of another.
    axon::Any
    soma::Any

    end

    """

    the meaning: i s the relationship between the concept

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
module ClimateERA

# Main file for the ClimateERA module that downloads and processes ECMWF
# reanalysis data.

## Modules Used
using Crayons
using Dates
using DelimitedFiles
using GeoRegions
using Glob
using JLD2
using Logging
using NCDatasets
using Printf
using Statistics

## Exporting the following functions:
export
        erastartup, erainitialize, eradownload, eraanalysis, eracompile, eralvlsplit,
        erasubregion, erancread,
        erarawname, erarawread, erarawfolder, erarawsave,
        eraananame, eraanaread, eraanafolder,
        eracmpname, eracmpread,
        erafolder, eraregfolder, eravarfolder, eraimgfolder,
        erarawregion, erarawpoint, erarawgrid,
        eraparameterload, eraparametercopy, eraparameteradd, erapressureload,
        eraregionstep,
        eplotsubregion, eplotlsm, eplotlsmgrid,
        queryeaction, queryedataset, queryemod, queryepar, putinfo,
        hrstep, hrindy

## Including other files in the module

include("startup.jl")
include("initialize.jl")
include("download.jl")
include("analysis.jl")
include("compile.jl")
#include("calculate.jl")
include("eraplots.jl")

include("frontend.jl")
include("backend.jl")
include("query.jl")

end # module