text:
# Global variables
const epsilon = 10.0^(-10);
const epsilon2 = 1e-8;
const dx=0.15
const dt = 0.04;

const total_rate = 40;
const maximum_stim=3 ### update this with the magnitude of the maximum stimulus ####
using Polynomials
using Distributions
using ForwardDiff
using Polynomials
import ForwardDiff.DiffBase

# 3rd party
# import Base.convert
# using MAT
# using ForwardDiff
# using Optim

# import ForwardDiff.DiffBase
# include("/home/genis/Desktop/parameter_functions.jl")
# include("/home/genis/model_fitting_julia/scripts/functions/utils.jl")

convert(::Type{Float64}, x::ForwardDiff.Dual) = Float64(x.value)
function convert(::Array{Float64}, x::Array{ForwardDiff.Dual})
    y = zeros(size(x));
    for i in 1:prod(size(x))
        y[i] = convert(Float64, x[i])
    end
    return y
end

immutable NumericPair{X,Y} <: Number
  x::X
  y::Y
end
Base.isless(a::NumericPair, b::NumericPair) = (a.x<b.x) || (a.x==b.x && a.y<b.y)

"""
function bin_centers = make_bins(B, dx, binN)
Makes a series of points that will indicate bin centers. The first and
last points will indicate sticky bins. No "bin edges" are made-- the edge
between two bins is always implicity at the halfway point between their
corresponding centers. The center bin is always at x=0; bin spacing
(except for last and first bins) is always dx; and the position
of the first and last bins is chosen so that |B| lies exactly at the
midpoint between 1st (sticky) and 2nd (first real) bins, as well as
exactly at the midpoint between last but one (last real) and last
(sticky) bins.
Playing nice with ForwardDiff means that the *number* of bins must be predetermined.
So this function will not actually set the number of bins; what it'll do is determine their
locations. To accomplish this separation, the function uses as a third parameter binN,
which should be equal to the number of bins with bin centers > 0, as follows:
   binN = ceil(B/dx)
and then the total number of bins will be 2*binN+1, with the center one always corresponding
to position zero. Use non-differentiable types for B and dx for this to work.
"""
function make_bins{T}(bins, B, dx::T, binN)
    cnt = 1
    for i=-binN:binN
        #println(cnt)
        bins[cnt] = i*dx
        cnt = cnt+1
    end

    if binN*dx == B
        bins[end] = B + dx
        bins[1] = -B - dx
    else
        bins[end] = 2*B - (binN-1)*dx
        bins[1] = -2*B + (binN-1)*dx
    end
end

function Fmatrix2{T}(F::AbstractArray{T,2},params::Vector, bin_centers,absorbing=0)

    stim=params[1]
    c2=params[2]
    c4=params[3]
    c6=params[4]
    sigma2=params[5]
    bias_d=params[6]
    hb_d=params[7]
    # println("sigmaaaa ",sigma2)
    # sigma2_sbin = convert(typeof(sigma2), sigma2)
    sigma2_sbin = convert(Float64, sigma2)

    # println(typeof(sigma2))
    # println(size(sigma2))
    # println(sigma2)
    # println("Converted ", sigma2_sbin)
    # println(typeof(sigma2_sbin))

    n_sbins = max(70, ceil(10*sqrt(sigma2_sbin)/dx))
    #println(n_sbins,sigma2_sbin)
    if absorbing==1
        F[1,1] = 1;
        F[end,end] = 1;
    else
        F[2,1] = 1;
        F[end-1,end] = 1;
    end
    swidth = 5.*sqrt(sigma2_sbin)

    sbinsize = swidth/n_sbins #sbins[2] - sbins[1]
    base_sbins    = collect(-swidth:sbinsize:swidth)

    ps       = exp.(-base_sbins.^2/(2*sigma2)) # exp(Array) -> exp.(x)
    ps       = ps/sum(ps);

    sbin_length = length(base_sbins)
    binN = length(bin_centers)

    mu = 0.
    for j in 2:binN-1
        mu = bin_centers[j] - (-(stim+bias_d+hb_d)-c2*bin_centers[j]+c4*bin_centers[j]^3+c6*bin_centers[j]^5)*dt#(exp(lam*dt))
        # mu=bin_centers[j] +stim*dt

        for k in 1:sbin_length
            sbin = (k-1)*sbinsize + mu - swidth

            if sbin <= bin_centers[1] #(bin_centers[1] + bin_centers[2])/2
                F[1,j] = F[1,j] + ps[k]
            elseif bin_centers[end] <= sbin#(bin_centers[end]+bin_centers[end-1])/2 <= sbins[k]
                F[end,j] = F[end,j] + ps[k]
            else # more condition
                if (sbin > bin_centers[1] && sbin < bin_centers[2])
                    lp = 1; hp = 2;
                elseif (sbin > bin_centers[end-1] && sbin < bin_centers[end])
                    lp = binN-1; hp = binN;
                else
                    lp = floor(Int,((sbin-bin_centers[2])/dx)) + 2#find(bin_centers .<= sbins[k])[end]
                    hp = ceil(Int,((sbin-bin_centers[2])/dx)) + 2#lp+1#Int(ceil((sbins[k]-bin_centers[2])/dx) + 1);
                end

                if lp == hp
                    F[lp,j] = F[lp,j] + ps[k]
                else
                    F[hp,j] = F[hp,j] + ps[k]*(sbin - bin_centers[lp])/(bin_centers[hp] - bin_centers[lp])
                    F[lp,j] = F[lp,j] + ps[k]*(bin_centers[hp] - sbin)/(bin_centers[hp] - bin_centers[lp])
                end
            end
        end
    end
end

function Fmatrix_simple{T}(F::AbstractArray{T,2},params::Vector, bin_centers,absorbing=0)

    stim=params[1]
    c2=params[2]
    c4=params[3]
    c6=params[4]
    sigma2=params[5]
    bias_d=params[6]
    hb_d=params[7]
    # println("sigmaaaa ",sigma2)
    # sigma2_sbin = convert(typeof(sigma2), sigma2)
    sigma2_sbin = convert(Float64, sigma2)

    # println(typeof(sigma2))
    # println(size(sigma2))
    # println(sigma2)
    # println("Converted ", sigma2_sbin)
    # println(typeof(sigma2_sbin))

    n_sbins = max(70, ceil(10*sqrt(sigma2_sbin)/dx))
    if absorbing==1
        F[1,1] = 1;
        F[end,end] = 1;
    else
        F[2,1] = 1;
        F[end-1,end] = 1;
    end
    swidth = 5.*sqrt(sigma2_sbin)

    sbinsize = swidth/n_sbins #sbins[2] - sbins[1]
    base_sbins    = collect(-swidth:sbinsize:swidth)

    ps       = exp.(-base_sbins.^2/(2*sigma2)) # exp(Array) -> exp.(x)
    ps       = ps/sum(ps);

    sbin_length = length(base_sbins)
    binN = length(bin_centers)

    mu = 0.
    for j in 2:binN-1

        #mu = bin_centers[j] - (-(stim+bias_d+hb_d)-c2*bin_centers[j]+c4*bin_centers[j]^3+c6*bin_centers[j]^5)*dt#(exp(lam*dt))
        aux_mu=stim+bias_d+hb_d
        mu=(bin_centers[j] + aux_mu/c2)*exp(c2*dt) - aux_mu/c2

        for k in 1:sbin_length
            sbin = (k-1)*sbinsize + mu - swidth

            if sbin <= bin_centers[1] #(bin_centers[1] + bin_centers[2])/2
                F[1,j] = F[1,j] + ps[k]
            elseif bin_centers[end] <= sbin#(bin_centers[end]+bin_centers[end-1])/2 <= sbins[k]
                F[end,j] = F[end,j] + ps[k]
            else # more condition
                control=true
                i=2
                while control && i< length(bin_centers)
                    if sbin< bin_centers[i]+dx/2
                        F[i,j]=F[i,j]+ps[k]
                        control=false
                    end
                    i=i+1
                end
            end
        end
    end
end

function ProbRight2_Crank(stim,history_bias,Tframe,param,model_f)

    Ntrials,Nframes=(size(stim))
    PR_model=zeros(Ntrials)
    Nsteps=convert(Int,floor(Tframe/dt))

    sigma_a,c2,c4,c6,bias_d,hbias_d,x0,ro,bins,Nbins,absorbing=param_aux(stim,param,model_f,Ntrials,Nframes,history_bias)
    s_trans=zeros(typeof(sigma_a),Nframes)

    for itrial in 1:Ntrials
        model_f[:st_f](stim[itrial,:],s_trans,param["st"])

        PR_model[itrial]=ProbRight_1stim_Crank_Nicolson(s_trans,history_bias[itrial],bins,Nbins,Nsteps,Nframes,absorbing,c2,c4,c6,bias_d,
        hbias_d,sigma_a,ro,x0[itrial,:])

    end
    return PR_model
end

function ProbRight_1stim_Crank_Nicolson(s_trans,hb,bins,Nbins,Nsteps,Nframes,absorbing,
    c2,c4,c6,bias_d, hbias_d,sigma_a,ro,x0)

    # P[:]=0
    # mu[:]=0
    mu=zeros(typeof(sigma_a),Nbins)
    P=zeros(typeof(sigma_a),Nframes*Nsteps+1,Nbins)
    P[1,:]=x0
    aux=dt/(2*dx)
    D=(sigma_a^2)/2

    nu=D*dt/(dx^2)

    M=Tridiagonal(zeros(Nbins-1),2*nu*ones(Nbins),zeros(Nbins-1))
    I=Diagonal(ones(Nbins))
    it_aux=1
    #println("D: ",D,"nu: ",nu)
    for iframe in 1:Nframes
        for ix in 1:Nbins
            #print(mu[ix]," ",s_trans[iframe]," ",bias_d[iframe]," ",hbias_d[iframe]," ",bins[ix]," ",c2[iframe]," ",c4[iframe]," ",c6[iframe],"\n")
            mu[ix]=-(-(s_trans[iframe]+bias_d[iframe]+hbias_d[iframe]*hb[1])-c2[iframe]*bins[ix]+c4[iframe]*bins[ix]^3+c6[iframe]*bins[ix]^5)

        end
        chi=(dt/dx)*mu
        M.dl[:]=-nu-0.5*chi[2:end]
        M.du[:]=-nu+0.5*chi[2:end]

        #build Matrix

        for it in 1:Nsteps
                #println(it_aux)
                P[it_aux+1,:]=(I+0.5*M) \ ( ( I-0.5*M)*P[it_aux,:]  )
                # indx=find(x->x<0,P[it_aux+1,:])
                # P[it_aux+1,indx]=0
                P[it_aux+1,:]=P[it_aux+1,:]/(sum(P[it_aux+1,:]))
                it_aux=it_aux+1

        end
    end

    pr=dot(P[end,:],ro)
    return pr
end

function ProbRight2_parallel(stim,history_bias,Tframe,param,model_f)

    Ntrials,Nframes=(size(stim))
    #PR_model=zeros(Ntrials)
    PR_model = SharedArray{ typeof(x[1]) }(Ntrials)

    Nsteps=convert(Int,floor(Tframe/dt))
    sigma_a,c2,c4,c6,bias_d,hbias_d,x0,ro,bins,Nbins,absorbing=param_aux(stim,param,model_f,Ntrials,Nframes,history_bias)
    println("Nbins",Nbins," bins[1] ",bins[1])

    @sync @parallel for itrial in 1:Ntrials
        F=zeros(typeof(sigma_a),Nbins,Nbins)
        s_trans=zeros(typeof(sigma_a),Nframes)

        model_f[:st_f](stim[itrial,:],s_trans,param["st"])

        PR_model[itrial]=ProbRight_1stim_fast(s_trans,history_bias[itrial],bins,Nbins,Nsteps,Nframes,absorbing,c2,c4,c6,bias_d,
        hbias_d,sigma_a,ro,x0[itrial,:],F)

    end
    return PR_model
end

function ProbRight2(stim,history_bias,Tframe,param,model_f)

    Ntrials,Nframes=(size(stim))
    PR_model=zeros(Ntrials)
    Nsteps=convert(Int,floor(Tframe/dt))
    sigma_a,c2,c4,c6,bias_d,hbias_d,x0,ro,bins,Nbins,absorbing=param_aux(stim,param,model_f,Ntrials,Nframes,history_bias)
    println("Nbins",Nbins," bins[1] ",bins[1])
    s_trans=zeros(typeof(sigma_a),Nframes)
    F=zeros(typeof(sigma_a),Nbins,Nbins)

    for itrial in 1:Ntrials
        model_f[:st_f](stim[itrial,:],s_trans,param["st"])

        PR_model[itrial]=ProbRight_1stim_fast(s_trans,history_bias[itrial],bins,Nbins,Nsteps,Nframes,absorbing,c2,c4,c6,bias_d,
        hbias_d,sigma_a,ro,x0[itrial,:],F)

    end
    return PR_model
end

function ProbRight_1stim_fast(s_trans,hb,bins,Nbins,Nsteps,Nframes,absorbing,
    c2,c4,c6,bias_d,hbias_d,sigma_a,ro,x0,F)

    F[:]=0
    #F=zeros(typeof(sigma_a),Nbins,Nbins)
    sigma_effective=dt*sigma_a^2   #here I transform std to variance. The energy landscape is symtric for pos and neg sitgma_a
    #x0=zeros(typeof(sigma_a),Nbins)
    #binN=convert(Int,(Nbins-1)/2)
    #x0[convert(Int,binN+1)]=1
    params2=zeros(typeof(c2[1]),7)
    params2[1]=s_trans[1]
    params2[2]=c2[1]
    params2[3]=c4[1]
    params2[4]=c6[1]
    params2[5]=sigma_effective
    params2[6]=bias_d[1]
    params2[7]=hbias_d[1]*hb
    # params2=[k*s[1],c2,c4,sigma_effective]
    #Fmatrix_simple(F,params2,bins,absorbing)
    Fmatrix2(F,params2,bins,absorbing)

    x=F*x0

    for iframe in 1:Nframes

        params2[1]=s_trans[iframe]
        params2[2]=c2[iframe]
        params2[3]=c4[iframe]
        params2[4]=c6[iframe]
        params2[6]=bias_d[iframe]
        params2[7]=hbias_d[iframe]*hb
        # params2=[k*s[iframe],c2,c4,sigma_effective]
        #F=zeros(typeof(sigma_a),Nbins,Nbins)

        if iframe==1
            N_aux=Nsteps-1
        else
            N_aux=Nsteps
            F[:]=0
            Fmatrix2(F,params2,bins,absorbing)
        end

        #println(N_aux)
        for iN in 1:N_aux
            x=F*x
        end
    end
    pr=dot(x,ro)

    return pr
end

# function ProbRight_1stim_fast_forward_euler(s_trans,hb,bins,Nbins,Nsteps,Nframes,absorbing;
#     c2=zeros(Nframes),c4=zeros(Nframes),c6=zeros(Nframes),bias_d=zeros(Nframes),
#     hbias_d=zeros(Nframes),sigma_a=0.,ro=vcat(zeros((Nbins-1)/2),[0.5],ones( (Nbins-1)/2)),
#     x0=vcat(zeros((Nbins-1)/2),[1.],zeros( (Nbins-1)/2))
#     )
#     Nframes=length(s_trans)
#     mu=zeros(typeof(sigma_a),Nbins)
#     P=zeros(typeof(sigma_a),Nbins,Nframes*Nsteps+1)
#
#     P[:,1]=x0
#     aux=dt/(2*dx)
#     D=(sigma_a^2)/2
#     aux2=D*dt/(dx^2)
#     for iframe in 1:Nframes
#         it_aux=(iframe-1)*Nsteps
#         for ix in 1:Nbins
#             #print(mu[ix]," ",s_trans[iframe]," ",bias_d[iframe]," ",hbias_d[iframe]," ",bins[ix]," ",c2[iframe]," ",c4[iframe]," ",c6[iframe],"\n")
#             mu[ix]=-(-(s_trans[iframe]+bias_d[iframe]+hbias_d[iframe]*hb[1])-c2[iframe]*bins[ix]+c4[iframe]*bins[ix]^3+c6[iframe]*bins[ix]^5)
#
#         end
#         for it in 1:Nsteps
#             for ix in 2:Nbins-1
#                 P[ix,it+1+it_aux]=P[ix,it+it_aux]-aux*( mu[ix+1]*P[ix+1,it+it_aux]-mu[ix-1]*P[ix-1,it+it_aux])+ aux2*(P[ix+1,it+it_aux]+P[ix-1,it+it_aux]-2*P[ix,it+it_aux])
#             end
#
#             ix=1
#             P[2,it+it_aux]=P[2,it+it_aux]-aux*( mu[ix+1]*P[ix+1,it+it_aux])+ aux2*(P[ix+1,it+it_aux]-2*P[ix,it+it_aux])
#             ix=Nbins-1
#             P[Nbins-2,it+it_aux]=P[Nbins-2,it+it_aux]-aux*(-mu[ix-1]*P[ix-1,it+it_aux])+ aux2*(P[ix-1,it+it_aux]-2*P[ix,it+it_aux])
#         end
#     end
#
#     pr=dot(P[:,end],ro)
#
#     return pr
# end

function LogLikelihood2(s,d,hb,bins,Nbins,Nsteps,Nframes,absorbing,model_f,
    c2,c4,c6,bias_d,hbias_d,sigma_a,ro,x0)

    F=zeros(typeof(c2[1]),Nbins,Nbins)
    pr=ProbRight_1stim_fast(s,hb,bins,Nbins,Nsteps,Nframes,absorbing,c2,c4,c6,bias_d,
    hbias_d,sigma_a,ro,x0,F)

    if d==1
        # return log(ProbRight_1stim(s,bins,Nsteps,Nframes; c2=c2,c4=c4,sigma_a=sigma_a,lapse=lapse,B=B,k=k)+epsilon)
        # pr=ProbRight_1stim(s,bins,Nsteps,Nframes; kwargs...)+epsilon2
        #println("pr3:", pr)
        return log( maximum( [pr epsilon2] ) )

    else
        #println("pr2:", pr)
        # pr=ProbRight_1stim(s,bins,Nsteps,Nframes; kwargs...)+epsilon2
        return log( maximum([1.0-pr epsilon2] ) )

        # return log(1-ProbRight_1stim(s,bins,Nsteps,Nframes; c2=c2,c4=c4,sigma_a=sigma_a,lapse=lapse,B=B,k=k)+epsilon)
    end
end

function LogLikelihood_sigma0(s,hb,bins,Nbins,Nsteps,Nframes,absorbing,model_f,
    c2,c4,c6,bias_d,hbias_d,sigma_a,ro,x0)
    F=zeros(typeof(c2[1]),Nbins,Nbins)

    pr=ProbRight_1stim_fast(s,hb,bins,Nbins,Nsteps,Nframes,absorbing,c2,c4,c6,bias_d,
    hbias_d,sigma_a,ro,x0,F)

    return log( maximum( [pr epsilon2] ) ), log( maximum([1.0-pr epsilon2] ) )

end

function param_aux(stim,param,model_f,Ntrials,Nframes,history_bias)

    sigma_a=param["sigma_a"][1]

    c2=zeros(typeof(param["st"][1]),Nframes) #We need to define here c2, array can not be define in functions because of ForwardDiff
    c4=zeros(typeof(param["st"][1]),Nframes)
    c6=zeros(typeof(param["st"][1]),Nframes)
    bias_d=zeros(typeof(param["st"][1]),Nframes)
    hbias_d=zeros(typeof(param["st"][1]),Nframes)
    model_f[:c2_f](c2,param["c2"])
    model_f[:c4_f](c4,param["c4"])
    if "c6" in keys(param)
        model_f[:c6_f](c6,param["c6"])
    else
        model_f[:c6_f](c6)
    end

    if "bias_d" in keys(param)
        model_f[:bias_d_f](bias_d,param["bias_d"])
    else
        model_f[:bias_d_f](bias_d)
    end

    if "hbias_d" in keys(param)
        model_f[:hbias_d_f](hbias_d,param["hbias_d"])
    else
        model_f[:hbias_d_f](hbias_d)
    end
    max_mu=zeros(1)
    stim_max=zeros(1)+maximum_stim
    model_f[:st_f](stim_max,max_mu,param["st"])

    if "B" in keys(param)
        B=param["B"]
        absorbing=param["absorbing"]
    else
        B,absorbing=compute_bound_all(max_mu[1],c2,c4,c6,sigma_a)
    end
    #println("B:",B[1])
    binN=Int(ceil(B[1]/dx))
    Nbins=2*binN+1
    bins=zeros(typeof(param["st"][1]),Nbins)
    make_bins(bins,B[1],dx,binN)

    ro=zeros(typeof(param["st"][1]),Nbins)
    x0=zeros(typeof(param["st"][1]),(Ntrials,Nbins))

    if model_f[:ro_f]==read_out_perfect
        model_f[:ro_f](ro,binN)
    elseif ro_f==read_out_TW
        model_f[:ro_f](ro,bins,c2[end],c4[end],c6[end])
    else:
        model_f[:ro_f](ro,bins,param["ro"])
    end

    if model_f[:x0_f]==initial_condition
        x0_a=zeros(typeof(param["st"][1]),Nbins)
        model_f[:x0_f](x0_a)
        x0[:]=repeat(x0_a,1,Ntrials)
    elseif model_f[:x0_f]==initial_condition_bias
        x0_a=zeros(typeof(param["st"][1]),Nbins)
        model_f[:x0_f](x0_a,bins,param["x0"])
        x0[:]=repeat(x0_a,1,Ntrials)
    elseif model_f[:x0_f]==initial_condition_bias_hbias
        x0_m=param["x0"][1]+param["x0"][2]*history_bias
        model_f[:x0_f](x0,bins,x0_m)
    end

    return sigma_a,c2,c4,c6,bias_d,hbias_d,x0,ro,bins,Nbins,absorbing

end

function ComputeLL2(stim,i_sigma,choices,history_bias,Tframe,args,x,consts,y,model_f;fastsigma0=true)

    Ntrials,Nframes=(size(stim))
    Nsteps=convert(Int,floor(Tframe/dt))
    param=make_dict(args,x,consts,y)

    sigma_a,c2,c4,c6,bias_d,hbias_d,x0,ro,bins,Nbins,absorbing=param_aux(stim,param,model_f,Ntrials,Nframes,history_bias)
    LLs = SharedArray{ typeof(x[1]) }(Ntrials)
    #LLs=zeros(typeof(sigma_a),Ntrials)
    s_trans=zeros(typeof(param["st"][1]),Nframes)

    ### comment if there is no sigma=0##
    if fastsigma0
        isig0=findfirst(x->x==0,i_sigma)
        s_r=abs.(stim[isig0,:])
        model_f[:st_f](s_r,s_trans,param["st"])
        ll_rr,ll_rl=LogLikelihood_sigma0(s_trans,0.0,bins,Nbins,Nsteps,Nframes,absorbing,model_f,
        c2,c4,c6,bias_d,hbias_d,sigma_a,ro,x0[1,:])
        #ll_rr,ll_rl ll of a right trial with sigma=0 and  decision r(l)
        s_l=-abs.(stim[isig0,:]) #left stimulus with sigma=0
        model_f[:st_f](s_l,s_trans,param["st"])
        ll_lr,ll_ll=LogLikelihood_sigma0(s_trans,0.0,bins,Nbins,Nsteps,Nframes,absorbing,model_f,
        c2,c4,c6,bias_d,hbias_d,sigma_a,ro,x0[1,:])
    end

    @sync @parallel for itrial in 1:Ntrials
    #for itrial in 1:Ntrials
        if i_sigma[itrial]==0
            if stim[itrial,1]>0
                if choices[itrial]==1
                    LLs[itrial]=ll_rr
                else
                    LLs[itrial]=ll_rl
                end
            else
                if choices[itrial]==1
                    LLs[itrial]=ll_lr
                else
                    LLs[itrial]=ll_ll
                end
            end

        else
            model_f[:st_f](stim[itrial,:],s_trans,param["st"])
            LLs[itrial]=LogLikelihood2(s_trans,choices[itrial],history_bias[itrial],bins,Nbins,Nsteps,Nframes,absorbing,model_f,
            c2,c4,c6,bias_d,hbias_d,sigma_a,ro,x0[itrial,:])

        end
    end

    LL = -sum(LLs)

    return  LL + model_f[:lambda]*sum(abs.(x))

end

function ComputeLL2_rats(stim,i_sigma,i_mu,mu_vec,choices,history_bias,Tframe,args,x,consts,y,model_f;)
    Ntrials,Nframes=(size(stim))
    Nsteps=convert(Int,floor(Tframe/dt))
    param=make_dict(args,x,consts,y)

    sigma_a,c2,c4,c6,bias_d,hbias_d,x0,ro,bins,Nbins,absorbing=param_aux(stim,param,model_f,Ntrials,Nframes,history_bias)
    LLs = SharedArray{ typeof(x[1]) }(Ntrials)
    #LLs=zeros(typeof(sigma_a),Ntrials)

    ll_r = zeros(typeof(sigma_a),size(mu_vec))
    ll_l = zeros(typeof(sigma_a),size(mu_vec))
    counter = 1
    s_trans=zeros(typeof(sigma_a),Nframes)
    aux = zeros(typeof(sigma_a),Nframes)
    for mu in mu_vec
        # i_sig0 = find(x->x==0, i_sigma)
        # x_sig0 = x0[i_sig0,:]
        # i_mu0 = i_mu[i_sig0]
        # imu=findfirst(x->x==mu,i_mu0)
        # if imu == 0
        #     # println("imu = 0")
        #     ll_r[counter],ll_l[counter] = 0,0
        # else
        # println("imu != 0")
        # x0_0 = x_sig0[imu,:] # mirar si sempre és la mateixa
        s = aux+mu
        # println(s)

        model_f[:st_f](s,s_trans,param["st"])
        ll_r[counter],ll_l[counter]=LogLikelihood_sigma0(s_trans,0.0,bins,Nbins,Nsteps,Nframes,absorbing,model_f,
        c2,c4,c6,bias_d,hbias_d,sigma_a,ro,x0[1,:])
        counter += 1
    end
    # s_trans=zeros(typeof(sigma_a),Nframes)
    @sync @parallel for itrial in 1:Ntrials
        # print(itrial)
    # for itrial in 1:Ntrials
        if i_sigma[itrial]==0
            # println("isigma = 0")
            mu_trial = i_mu[itrial]
            mu_ind = findfirst(x -> x == mu_trial,mu_vec) # Provar de fer entrar un vector d'index mu
            if choices[itrial] == 1
                LLs[itrial] = ll_r[mu_ind[1]]
            else
                LLs[itrial] = ll_l[mu_ind[1]]
            end

        else
            # println("isigma != 0")
            model_f[:st_f](stim[itrial,:],s_trans,param["st"])
            LLs[itrial]=LogLikelihood2(s_trans,choices[itrial],history_bias[itrial],bins,Nbins,Nsteps,Nframes,absorbing,model_f,
            c2,c4,c6,bias_d,hbias_d,sigma_a,ro,x0[itrial,:])
        end
    end

    LL = -sum(LLs)

    return  LL +model_f[:lambda]*sum(abs.(x))

end

function compute_bound_all(max_mu,c2,c4,c6,sigmas)
    y=zeros(Float64,4)
    b=zeros(Float64,length(c2))
    sign=zeros(Float64,length(c2))
    for iframe in 1:length(c2)
        y[1]=c2[iframe]
        y[2]=c4[iframe]
        y[3]=c6[iframe]
        #y[4]=sigmas[iframe]
        y[4]=sigmas[1]

        b[iframe],sign[iframe]=compute_bound(;max_mu=max_mu,c2=y[1],c4=y[2],c6=y[3],sigma_a=y[4])
    end
    idx=find(b.==maximum(b))[1]
    return b[idx],sign[idx]
end

function compute_bound(;max_mu=0,c2=1,c4=1,c6=0,sigma_a=0.1)
    # if sigma_a<2.
    #     a=Polynomials.poly([,1.035,5.471])
    # else
    #     a=Polynomials.poly([34.626,-23.743,10.040])
    # end
    #a=Polynomials.Poly([0.0038,2.575,1.30375])### 99%
    a=Polynomials.Poly([0.6865, 2.6837, 2.6068])### 99,9%
    fprima=Polynomials.polyval(a,sigma_a)
    #p=Poly([1.64*2./sqrt(dt)*sigma_a^2, -c2, 0, c4,0,c6])
    p=Poly([-fprima-max_mu,-c2,0,c4,0,c6])
    roots_p=roots(p)
    # println("roots :",roots_p,c2," ", c4, " ",sigma_a)
    max_bound=40.
    if length(roots_p)==0
        b=max_bound

    elseif typeof(roots_p[1])!=Float64
        aux=filter(x-> x.im==0.0,roots_p)
        roots_real=zeros(length(aux))
        for iroot in 1:length(aux)
            roots_real[iroot]=abs(aux[iroot].re)
        end
        b=minimum( [maximum(roots_real),max_bound])
    else
        b=minimum( [maximum(abs.(roots_p)),max_bound])
    end
    # println(roots_p," ",b)
    return b+dx,-sign(b*potential_derivative(b,0.,c2,c4,c6))
end

function potential_derivative(x,s,c2,c4,c6)

    return -s-c2*x+c4*x^3+c6*x^5
end

function ModelFitting(stim,i_sigma,choices,history_bias, Tframe::Float64, args,consts,ub,lb, x_ini::Vector{Float64},y,model_f;fastsigma0=true)

    function LL_f(x)
        # println("param: ",x)ComputeLL2
        return ComputeLL2(stim,i_sigma,choices,history_bias,Tframe, args,x,consts,y,model_f,fastsigma0=fastsigma0)
    end

    println("x_ini hola ",x_ini)
    #println("lower bound hola ",lb)
    #println("upper bound hola ",ub)
    d4 = OnceDifferentiable(LL_f,x_ini;autodiff=:forward)
    # d4 = OnceDifferentiable(Optim.only_fg!(LL_fg!),x_ini)
    #
    # println(d4.f(x_ini))
    # grads=zeros(4)
    # # d4.df(x_ini,grads)
    # d4.df(grads,x_ini)
    # println(grads)

    tic()

    history = optimize(d4, lb, ub, x_ini, Fminbox(LBFGS()),Optim.Options(show_trace = false,store_trace=false,extended_trace = false,iterations=50,g_tol = 1e-8,f_tol = 1e-8,x_tol=1e-8) )
    # history = optimize(d4, x_ini, LBFGS(),Optim.Options(store_trace=true,extended_trace = true,time_limit=2*3600,show_trace = true,iterations=50,g_tol = 1e-5,f_tol = 1e-4,x_tol=1e-5,g_calls_limit=2000) )
    # history = optimize(d4, x_ini, LBFGS(),Optim.Options(store_trace=true,extended_trace = true,time_limit=2*3600,show_trace = true,iterations=50,g_tol = 1e-5,f_tol = 1e-4,x_tol=1e-5,g_calls_limit=2000) )

    # history = optimize(d4, l, u, x_ini, Fminbox(LBFGS()),Optim.Options(g_tol = 1e-5,
    #                                                                        iterations=1,
    #                                                                        x_tol = 1e-2,
    #                                                                        f_tol = 1e-1,                                                                        iterations = 10,
    #                                                                        store_trace = true,
    #                                                                        show_trace = true,
    #                                                                        extended_trace = true,
    #                                                                        ))

    # history = optimize(d4, l, u, x_ini, Fminbox(LBFGS()),Optim.Options(g_tol = 1e-5,
    #                                                                        x_tol = 1e-5,
    #                                                                        f_tol = 1e-2,                                                                        iterations = 10,
    #                                                                        store_trace = true,
    #                                                                        show_trace = true,
    #                                                                        extended_trace = true,
    #                                                                        iterations=500))

     # ;
     #         optimizer_o = Optim.Options(g_tol = 1e-12,
     #                                                                        x_tol = 1e-10,
     #                                                                        f_tol = 1e-6,                                                                        iterations = 10,
     #                                                                        store_trace = true,
     #                                                                        show_trace = true,
     #                                                                        extended_trace = true))
    # history=optimize(LL_f,x_ini,LBFGS(); autodiff=:forward)

    # history=Optim.optimize(Optim.only_fg!(LL_fg!), x_ini, Optim.LBFGS())
    # history=Optim.optimize(Optim.only_fg!(LL_fg!), l, u, Fminbox(),x_ini, Optim.LBFGS())

    fit_time=toc()
    # println(history.minimizer)
    println(history)

    x_bf = history.minimizer

    Hess=ForwardDiff.hessian(LL_f, x_bf)
    Gs = zeros(length(history.trace),length(x_ini))
    Xs = zeros(length(history.trace),length(x_ini))
    fs = zeros(length(history.trace))

    ci=compute_confidence_interval(Hess)

    for i=1:length(history.trace)
        tt = getfield(history.trace[i],:metadata)
        fs[i] = getfield(history.trace[i],:value)
        Gs[i,:] = tt["g(x)"]
        Xs[i,:] = tt["x"]
    end
    param3=make_dict(args,history.minimizer,consts,y)
    D = Dict([("x_ini",x_ini),
                ("parameters",param3),
                ("trials",size(stim)[1]),
                ("f",history.minimum),
                ("x_converged",history.x_converged),
                ("f_converged",history.f_converged),
                ("g_converged",history.g_converged),
                ("grad_trace",Gs),
                ("f_trace",fs),
                ("x_trace",Xs),
                ("fit_time",fit_time),
                ("x_bf",history.minimizer),
                ("myfval", history.minimum),
                ("f_calls", history.f_calls),
                ("g_calls",history.g_calls),
                ("iterations",history.iterations),
                ("Hessian",Hess),
                ("ci",ci),
                ])
    return D

end

function compute_confidence_interval(H)
    z_aux=1.96
    ci=zeros(size(H)[1])
    try
        HI=inv(H)
        if all(diag(HI).>0)
            ci[:]=z_aux.*sqrt.(diag(HI))
            #println("A minimum ", ci)
        else
            ci[:]=-1.0
            #println("Not a minimum ", diag(HI))
        end

    catch e
        ci[:]=-2.0
        #println("H is not invertible")
    end
    return ci
end

function ModelFitting_rats(stim,i_sigma,i_mu,mu_vec,choices,history_bias, Tframe::Float64, args,consts,ub,lb, x_ini::Vector{Float64},y,model_f)

    function LL_f(x)
        # println("param: ",x)ComputeLL2
        return ComputeLL2_rats(stim,i_sigma,i_mu,mu_vec,choices,history_bias,Tframe, args,x,consts,y,model_f)
    end

    println("x_ini hola ",x_ini)
    #println("lower bound hola ",lb)
    #println("upper bound hola ",ub)
    d4 = OnceDifferentiable(LL_f,x_ini;autodiff=:forward)
    # d4 = OnceDifferentiable(Optim.only_fg!(LL_fg!),x_ini)
    #
    # println(d4.f(x_ini))
    # grads=zeros(4)
    # # d4.df(x_ini,grads)
    # d4.df(grads,x_ini)
    # println(grads)

    tic()

    history = optimize(d4, lb, ub, x_ini, Fminbox(LBFGS()),Optim.Options(show_trace = false,store_trace=false,extended_trace = false,iterations=500,g_tol = 1e-8,f_tol = 1e-8,x_tol=1e-8) )
    # history = optimize(d4, x_ini, LBFGS(),Optim.Options(store_trace=true,extended_trace = true,time_limit=2*3600,show_trace = true,iterations=50,g_tol = 1e-5,f_tol = 1e-4,x_tol=1e-5,g_calls_limit=2000) )
    # history = optimize(d4, x_ini, LBFGS(),Optim.Options(store_trace=true,extended_trace = true,time_limit=2*3600,show_trace = true,iterations=50,g_tol = 1e-5,f_tol = 1e-4,x_tol=1e-5,g_calls_limit=2000) )

    # history = optimize(d4, l, u, x_ini, Fminbox(LBFGS()),Optim.Options(g_tol = 1e-5,
    #                                                                        iterations=1,
    #                                                                        x_tol = 1e-2,
    #                                                                        f_tol = 1e-1,                                                                        iterations = 10,
    #                                                                        store_trace = true,
    #                                                                        show_trace = true,
    #                                                                        extended_trace = true,
    #                                                                        ))

    # history = optimize(d4, l, u, x_ini, Fminbox(LBFGS()),Optim.Options(g_tol = 1e-5,
    #                                                                        x_tol = 1e-5,
    #                                                                        f_tol = 1e-2,                                                                        iterations = 10,
    #                                                                        store_trace = true,
    #                                                                        show_trace = true,
    #                                                                        extended_trace = true,
    #                                                                        iterations=500))

     # ;
     #         optimizer_o = Optim.Options(g_tol = 1e-12,
     #                                                                        x_tol = 1e-10,
     #                                                                        f_tol = 1e-6,                                                                        iterations = 10,
     #                                                                        store_trace = true,
     #                                                                        show_trace = true,
     #                                                                        extended_trace = true))
    # history=optimize(LL_f,x_ini,LBFGS(); autodiff=:forward)

    # history=Optim.optimize(Optim.only_fg!(LL_fg!), x_ini, Optim.LBFGS())
    # history=Optim.optimize(Optim.only_fg!(LL_fg!), l, u, Fminbox(),x_ini, Optim.LBFGS())

    fit_time=toc()
    # println(history.minimizer)
    println(history)

    x_bf = history.minimizer

    Hess=ForwardDiff.hessian(LL_f, x_bf)
    Gs = zeros(length(history.trace),length(x_ini))
    Xs = zeros(length(history.trace),length(x_ini))
    fs = zeros(length(history.trace))

    for i=1:length(history.trace)
        tt = getfield(history.trace[i],:metadata)
        fs[i] = getfield(history.trace[i],:value)
        Gs[i,:] = tt["g(x)"]
        Xs[i,:] = tt["x"]
    end
    param3=make_dict(args,history.minimizer,consts,y)
    D = Dict([("x_ini",x_ini),
                ("parameters",param3),
                ("trials",size(stim)[1]),
                ("f",history.minimum),
                ("x_converged",history.x_converged),
                ("f_converged",history.f_converged),
                ("g_converged",history.g_converged),
                ("grad_trace",Gs),
                ("f_trace",fs),
                ("x_trace",Xs),
                ("fit_time",fit_time),
                ("x_bf",history.minimizer),
                ("myfval", history.minimum),
                ("f_calls", history.f_calls),
                ("g_calls",history.g_calls),
                ("iterations",history.iterations),
                ("Hessian",Hess),
                ])
    return D

end

function LL_function_grads(stim,i_sigma,choices,history_bias,Tframe, args, x,consts,y,model_f)

    function LL_f23(z)
        # println("param: ",x)
        return ComputeLL2(stim,i_sigma,choices,history_bias,Tframe, args,z,consts,y,model_f)
    end

    d4=OnceDifferentiable(LL_f23,x;autodiff=:forward)
    # @time f=d4.f(x_ini)
    grads=zeros( length(x))

    f=d4.fdf(grads,x)

    # d4.df(x_ini,grads)
    # @time d4.df(grads,x_ini)
    println("grad: ",grads)
    # tic()
    # Hess=ForwardDiff.hessian(LL_f23, x)
    # toc()
    return f,grads#,Hess

end

function LL_from_PR(PR,choices)
    LL=0.
    for i in 1:length(choices)
        if choices[i]==1
            LL=LL-log(PR[i]+epsilon2)
        else
            try
                LL=LL-log(1-PR[i]+epsilon2)
            catch e
                LL=LL
            end
        end
    end
    return LL

end

function pdf_x_time_1stim_forward_euler(s,hb,Tframe,param,model_f)

    Nframes=length(s)
    stim_trans=zeros(Nframes)
    # println(param)
    sigma_a=param["sigma_a"][1]
    c2=zeros(typeof(sigma_a),Nframes) #We need to define hwere c2, array can not be define in functions because of ForwardDiff
    c4=zeros(typeof(sigma_a),Nframes)
    c6=zeros(typeof(sigma_a),Nframes)
    bias_d=zeros(typeof(sigma_a),Nframes)
    hbias_d=zeros(typeof(sigma_a),Nframes)
    model_f[:c2_f](c2,param["c2"])
    model_f[:c4_f](c4,param["c4"])

    if "c6" in keys(param)
        model_f[:c6_f](c6,param["c6"])
    else
        model_f[:c6_f](c6)
    end

    if "bias_d" in keys(param)
        model_f[:bias_d_f](bias_d,param["bias_d"])
    else
        model_f[:bias_d_f](bias_d)
    end

    if "hbias_d" in keys(param)
        #println(hbias_d)
        model_f[:hbias_d_f](hbias_d,param["hbias_d"])
    else
        model_f[:hbias_d_f](hbias_d)
    end
    Nsteps=convert(Int,floor(Tframe/dt))
    #println("Nsteps ",Nsteps,"Nframes: ",Nframes)

    # y=zeros(Float64,4) ### ForwardDiff is blind to thwe computation of the bound. It is not important for LL computation ####################
    # y[1]=minimum(c2)
    # y[2]=c4[1] ###### This is assuming c4 is constant during stimulus ############
    # y[3]=c6[1]
    # y[4]=sigma_a
    # B,absorbing=compute_bound(c2=y[1],c4=y[2],c6=y[3],sigma_a=y[4])

    s_trans=zeros(Nframes)
    model_f[:st_f](s,s_trans,param["st"])
    max_mu=maximum(abs(s))

    B,absorbing=compute_bound_all(max_mu,c2,c4,c6,sigma_a)

    binN=Int(ceil(B/dx))
    Nbins=2*binN+1
    bins=zeros(Float64,Nbins)
    make_bins(bins,B,dx,binN)

    sigma_effective=dt*sigma_a^2   #here I transform std to variance. The energy landscape is symtric for pos and neg sitgma_a
    binN=convert(Int,(Nbins-1)/2)

    Ntrials=1
    x0_a=zeros(typeof(sigma_a),Nbins)
    if model_f[:x0_f]==initial_condition
        model_f[:x0_f](x0_a)
        x0=repeat(x0_a,1,Ntrials)
    elseif model_f[:x0_f]==initial_condition_bias
        x0_a=zeros(typeof(sigma_a),Nbins)
        model_f[:x0_f](x0_a,bins,param["x0"])
        x0=repeat(x0_a,1,Ntrials)
    elseif model_f[:x0_f]==initial_condition_bias_hbias
        x0=zeros(typeof(sigma_a),(Ntrials,Nbins))
        x0_m=param["x0"][1]+param["x0"][2]*hb
        model_f[:x0_f](x0,bins,x0_m)
    end

    Nx=Nbins
    mu=zeros(typeof(sigma_a),Nx)
    P=zeros(typeof(sigma_a),Nx,Nframes*Nsteps+1)

    P[:,1]=x0
    aux=dt/(2*dx)
    D=(sigma_a^2)/2
    aux2=D*dt/(dx^2)
    for iframe in 1:Nframes
        it_aux=(iframe-1)*Nsteps
        for ix in 1:Nx
            #print(mu[ix]," ",s_trans[iframe]," ",bias_d[iframe]," ",hbias_d[iframe]," ",bins[ix]," ",c2[iframe]," ",c4[iframe]," ",c6[iframe],"\n")
            mu[ix]=-(-(s_trans[iframe]+bias_d[iframe]+hbias_d[iframe]*hb[1])-c2[iframe]*bins[ix]+c4[iframe]*bins[ix]^3+c6[iframe]*bins[ix]^5)

        end
        for it in 1:Nsteps
            for ix in 2:Nx-1
                P[ix,it+1+it_aux]=P[ix,it+it_aux]-aux*( mu[ix+1]*P[ix+1,it+it_aux]-mu[ix-1]*P[ix-1,it+it_aux])+ aux2*(P[ix+1,it+it_aux]+P[ix-1,it+it_aux]-2*P[ix,it+it_aux])
            end

            ix=1
            P[2,it+it_aux]=P[2,it+it_aux]-aux*( mu[ix+1]*P[ix+1,it+it_aux])+ aux2*(P[ix+1,it+it_aux]-2*P[ix,it+it_aux])
            ix=Nx-1
            P[Nx-2,it+it_aux]=P[Nx-2,it+it_aux]-aux*(-mu[ix-1]*P[ix-1,it+it_aux])+ aux2*(P[ix-1,it+it_aux]-2*P[ix,it+it_aux])
        end
    end
    ro=vcat(zeros(binN),[0.5],ones(binN))
    print("PR: ", dot(P[:,end],ro) )

    return P,bins
end

function pdf_x_time_1stim(s,Tframe,param,model_f,history_bias;B=0)

    # Ntrials,Nframes=(size(s))
    Ntrials=1
    Nframes=size(s)[1]
    Nsteps=convert(Int,floor(Tframe/dt))
    #param=make_dict(args,x,consts,y)

    sigma_a,c2,c4,c6,bias_d,hbias_d,x0,ro,bins,Nbins,absorbing=param_aux(s,param,model_f,Ntrials,Nframes,history_bias)

    #
    # Nframes=length(s)
    # stim_trans=zeros(Nframes)
    # sigma_a=param["sigma_a"][1]
    # c2=zeros(typeof(sigma_a),Nframes) #We need to define hwere c2, array can not be define in functions because of ForwardDiff
    # c4=zeros(typeof(sigma_a),Nframes)
    # c6=zeros(typeof(sigma_a),Nframes)
    # model_f[:c2_f](c2,param["c2"])
    # model_f[:c4_f](c4,param["c4"])
    #
    # bias_d=zeros(typeof(sigma_a),Nframes)
    # hbias_d=zeros(typeof(sigma_a),Nframes)
    #
    # if "bias_d" in keys(param)
    #     model_f[:bias_d_f](bias_d,param["bias_d"])
    # else
    #     model_f[:bias_d_f](bias_d)
    # end
    #
    # if "hbias_d" in keys(param)
    #     model_f[:hbias_d_f](hbias_d,param["hbias_d"])
    # else
    #     model_f[:hbias_d_f](hbias_d)
    # end
    #
    #
    # if "c6" in keys(param)
    #     model_f[:c6_f](c6,param["c6"])
    # else
    #     model_f[:c6_f](c6)
    # end
    #
    #
    # Nsteps=convert(Int,floor(Tframe/dt))
    #
    #
    # if B==0
    #     max_mu=zeros(1)
    #     stim_max=zeros(1)+maximum_stim
    #     model_f[:st_f](stim_max,max_mu,param["st"])
    #
    #     B,absorbing=compute_bound_all(max_mu[1],c2,c4,c6,sigma_a)
    #     #compute_bound_all(c2,c4,c6,sigma_a)
    # else
    #     absorbing=1
    # end
    #

    s_trans=zeros(Nframes)
    model_f[:st_f](s,s_trans,param["st"])
    println("s_trans: ",s_trans)

    # binN=Int(ceil(B/dx))
    # Nbins=2*binN+1
    # bins=zeros(Float64,Nbins)
    # make_bins(bins,B,dx,binN)
    #

    #
    #
    sigma_effective=dt*sigma_a^2   #here I transform std to variance. The energy landscape is symtric for pos and neg sitgma_a
    # binN=convert(Int,(Nbins-1)/2)
    # Ntrials=1
    # itrial=1
    # x0_a=zeros(typeof(sigma_a),Nbins)
    # if model_f[:x0_f]==initial_condition
    #     model_f[:x0_f](x0_a)
    #     x0=repeat(x0_a,1,Ntrials)
    # elseif model_f[:x0_f]==initial_condition_bias
    #     x0_a=zeros(typeof(sigma_a),Nbins)
    #     model_f[:x0_f](x0_a,bins,param["x0"])
    #     x0=repeat(x0_a,1,Ntrials)
    # elseif model_f[:x0_f]==initial_condition_bias_hbias
    #     x0=zeros(typeof(sigma_a),(Ntrials,Nbins))
    #     x0_m=param["x0"][1]+param["x0"][2]*hb
    #     model_f[:x0_f](x0,bins,x0_m)
    # end

    #println("sigma",sigma_a," c2 ",c2,"c4 ",c4," bias_d ",bias_d)
    #println(" hbias_d ",hbias_d," bins ", bins," Nbins ", Nbins, "abs", absorbing)

    F=zeros(typeof(sigma_a),Nbins,Nbins)

    params2=zeros(typeof(sigma_a),7)
    params2[1]=s_trans[1]
    params2[2]=c2[1]
    params2[3]=c4[1]
    params2[4]=c6[1]
    params2[5]=sigma_effective
    params2[6]=bias_d[1]
    params2[7]=hbias_d[1]*history_bias[1]

    x=zeros(typeof(sigma_a),Nsteps*Nframes,Nbins)
    j=1
    itrial=1
    Fmatrix2(F,params2,bins,absorbing)
    x[1,:]=F*x0[itrial,:]

    for iframe in 1:Nframes

        params2[1]=s_trans[iframe]
        params2[2]=c2[iframe]
        params2[3]=c4[iframe]
        params2[4]=c6[iframe]
        F=zeros(typeof(sigma_a),Nbins,Nbins)
        Fmatrix2(F,params2,bins,absorbing)
        if iframe==1
            N_aux=Nsteps-1
        else
            N_aux=Nsteps
        end
        for iN in 1:N_aux
            #println(sum(x[j,:]))
            x[j+1,:]=F*x[j,:]
            j+=1
        end

    end
    return x,bins
end

######################## functions PI #################

function ProbRightPI(stim,history_bias,Tframe,param,model_f)

    Ntrials,Nframes=(size(stim))
    PR_model=zeros(typeof(param["st"][1]),Ntrials)
    T=Nframes*Tframe
    s_trans=zeros(typeof(param["st"][1]),Nframes)
    std=sqrt(T)*param["sigma_a"][1]
    for itrial in 1:Ntrials
        model_f[:st_f](stim[itrial,:],s_trans,param["st"])
        mean=param["x0"][1]+sum(Tframe*s_trans) +T*param["bias_d"][1]
        f=Normal(mean,std)
        PR_model[itrial]=1-cdf(f,0)
    end

    return PR_model
end

function ComputeLLPI(stim,choices,history_bias,Tframe,args,consts, x,y,model_f)

    param=make_dict(args,x,consts,y)
    PR=ProbRightPI(stim,history_bias,Tframe,param,model_f)
    LL=LL_from_PR(PR,choices)

    return LL+model_f[:lambda]*sum(abs.(x))

end

function ModelFittingPI(stim,choices,history_bias, Tframe, args,consts,ub,lb, x_ini,y,model_f)

    function LL_f(x)
        # println("param: ",x)ComputeLL2
        return ComputeLLPI(stim,choices,history_bias,Tframe, args,consts,x,y,model_f)
    end

    println("x_ini hola ",x_ini)

    d4 = OnceDifferentiable(LL_f,x_ini;autodiff=:forward)

    tic()
    history = optimize(d4, lb, ub, x_ini, Fminbox(LBFGS()),Optim.Options(show_trace = false,store_trace=false,extended_trace = false,iterations=500,g_tol = 1e-8,f_tol = 1e-8,x_tol=1e-8) )
    fit_time=toc()
    # println(history.minimizer)
    println(history)

    x_bf = history.minimizer

    Hess=ForwardDiff.hessian(LL_f, x_bf)
    Gs = zeros(length(history.trace),length(x_ini))
    Xs = zeros(length(history.trace),length(x_ini))
    fs = zeros(length(history.trace))

    ci=compute_confidence_interval(Hess)

    for i=1:length(history.trace)
        tt = getfield(history.trace[i],:metadata)
        fs[i] = getfield(history.trace[i],:value)
        Gs[i,:] = tt["g(x)"]
        Xs[i,:] = tt["x"]
    end
    param3=make_dict(args,history.minimizer,consts,y)
    D = Dict([("x_ini",x_ini),
                ("parameters",param3),
                ("trials",size(stim)[1]),
                ("f",history.minimum),
                ("x_converged",history.x_converged),
                ("f_converged",history.f_converged),
                ("g_converged",history.g_converged),
                ("grad_trace",Gs),
                ("f_trace",fs),
                ("x_trace",Xs),
                ("fit_time",fit_time),
                ("x_bf",history.minimizer),
                ("myfval", history.minimum),
                ("f_calls", history.f_calls),
                ("g_calls",history.g_calls),
                ("iterations",history.iterations),
                ("Hessian",Hess),
                ("ci",ci),
                ])
    return D

end

function LL_function_grads_PI(stim,choices,history_bias,Tframe, args,consts,x,y,model_f)

    function LL_f23(z)
        return ComputeLLPI(stim,choices,history_bias,Tframe, args,consts,z,y,model_f)
    end

    d4=OnceDifferentiable(LL_f23,x;autodiff=:forward)
    # @time f=d4.f(x_ini)
    grads=zeros( length(x))

    f=d4.fdf(grads,x)

    # d4.df(x_ini,grads)
    # @time d4.df(grads,x_ini)
    println("grad: ",grads)
    # tic()
    # Hess=ForwardDiff.hessian(LL_f23, x)
    # toc()
    return f,grads#,Hess

end

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
using ProgressMeter
using Printf
using MATLAB
using BSON: @save, @load
using Dates
include("Dimer_Lagrangian_Descriptor.jl")

@everywhere begin
    include("Dimer_Lagrangian_Descriptor.jl")
    push!(LOAD_PATH, "/home/brandon_behring/Desktop/Leap_Frog_Code")
    using ProgressMeter
    using Printf
    using MATLAB
    using .Dimer_Lagrangian_Descriptor
    using BSON: @save, @load
    using Dates
    # H=range(.11,stop=0.145,length=107)
    # Energy=0.112012
    # Energy=.125789
    # Energy=.12501
    # Energy=.12501
    # Energy=.1249
    # Energy=.125

    N=Int(1e4)
    Energy=.129

    Q_start=.1
    Q_end=.5
    count = 1
    dp=1.0e-7
    Q=range(Q_start, stop=Q_end, length=N)
    P=range(-dp,stop=dp, length=5)
    P=0.0
    Nx=length(P)
    mesh=[(0.0,q) for q in Q]
    mesh = [(p, q) for q in Q, p in P]
    location_fig="/media/brandon_behring/Extra_Space/MATLAB_FIGURES/"
    right_now = replace(replace(replace(string(Dates.now()),"."=>"_"),":"=>"_"),"-"=>"_")
    t=1:50
end
while count <= length(t)
    @everywhere t_end=t[count]
    @everywhere h=replace(@sprintf("%.15f",Energy),"."=>"_")
    @everywhere Ns = @sprintf("_%5d",N)
    @everywhere tend= @sprintf("_%d",t_end)
    println(tend)
    @everywhere file_name=location_fig* "ld_grad_test_"*h*Ns*tend*right_now*".fig"
    @everywhere tend_T= @sprintf("%d",t_end)
    @everywhere h_title =h*" with tend="*tend_T*" and "*Ns
LD=@showprogress pmap(Dimer_Lagrangian_Descriptor.Dimer_Lagrangian_Descriptor_Function, mesh, Energy*ones(N,Nx), t_end*ones(N,Nx) )
dx=(Q_end-Q_start)/N
dy=2*dp/5
# gradM=Dimer_Lagrangian_Descriptor.gradient_matrix_4(LD,dy,dx)
# gradM2=Dimer_Lagrangian_Descriptor.D_X( LD[:,3],dx)
gradM2=Dimer_Lagrangian_Descriptor.D_X( LD,dx)
# gradMl=zeros(length(gradM))
gradMx=zeros(length(gradM2))
for i=1:length(gradM2)
    # gradMl[i]=log10(gradM[i])
    gradMx[i]=log10(gradM2[i])
end

# d2M=log10.(Dimer_Lagrangian_Descriptor.D_X(gradM,dx))
Q1=Q[3:N-2]
LDn=5*( LD/maximum(filter(!isnan,LD))-.7*ones(length(LD)) )
mat"figure();set(gcf, 'Position',  [0, 0, 2500, 1500]); hold on;"
mat"title($h_title)"
mat"plot($Q, $LDn,'b')"
# mat"plot($Q1, $gradMl,'r')"
mat"plot($Q1, $gradMx,'r')"
mat"savefig($file_name)"

@everywhere global count += 1

end

ext: jl


text:
module ChebyshevPicardIteration

using LinearAlgebra
using StaticArrays

include("ChebyshevPolynomials.jl")
include("Integrator.jl")

export Integrate

end


text:
baseprice = (3,1,2,5)
currentprice = (5,3,1,6)
baseqty = (25,50,30,15)
currentqty = (28,60,30,12)

@test FisherIndex(baseprice, currentprice, baseqty, currentqty)==1.55835405667883


text:
# (-3, i) = long field
# (-2, i) = transverse field off-diag
# (-1, i) = transverse field diag
# (0, 0) = id
# (i, j) = diag bond op

function make_prob_vector(J::AbstractMatrix{T}, hx::AbstractVector{T}) where T
    ops = Vector{NTuple{2, Int}}(undef, 0)
    p = Vector{T}(undef, 0)

    k = 0
    for i in eachindex(hx)
        if hx[i] != 0
            push!(ops, (-1, i))
            push!(p, h)
        end
    end

    # only take J_ij terms from upper diagonal since it's symmetric
    for j in axes(J, 2), i in axes(J, 1)
        if i < j
            if J[i, j] != 0
                push!(ops, op)
                push!(p, 2*abs(J))
            end
        end
    end

    return ops, p
end

function make_prob_vector(bond_spins::Vector{NTuple{2,Int}}, Ns::Int, J::T, h::T) where T
    ops = Vector{NTuple{2, Int}}(undef, 0)
    p = Vector{T}(undef, 0)

    if !iszero(h)
        for i in 1:Ns
            push!(ops, (-1, i))
            push!(p, h)
        end
    end

    if !iszero(J)
        for op in bond_spins
            push!(ops, op)
            push!(p, 2*abs(J))
        end
    end

    return ops, p
end

###################

abstract type AbstractOperatorSampler{N, T, P <: AbstractProbabilityVector{T}} end
firstindex(::AbstractOperatorSampler) = 1
lastindex(os::AbstractOperatorSampler) = length(os)

struct OperatorSampler{N, T, P} <: AbstractOperatorSampler{N, T, P}
    operators::Vector{NTuple{N, Int}}
    pvec::P
end

# samples operators using a heap
# based on a blog post by Tim Vieira
# https://timvieira.github.io/blog/post/2016/11/21/heaps-for-incremental-computation/
function OperatorSampler(operators::Vector{NTuple{N, Int}}, p::Vector{T}) where {N, T <: Real}
    @assert length(operators) == length(p) "Given vectors must have the same length!"
    pvec = probability_vector(p)
    return OperatorSampler{N, T, typeof(pvec)}(operators, pvec)
end
rand(os::OperatorSampler{N}) where N = @inbounds os.operators[rand(os.pvec)]
@inline length(os::OperatorSampler) = length(os.operators)

function cluster_probs_vec(operators::Vector{NTuple{N, Int}}, p::AbstractVector{T}) where {N, T <: Real}
    perm = sortperm(p)
    p = p[perm]
    operators = operators[perm]

    uniq_p = T[]
    uniq_ops = Vector{Vector{NTuple{N, Int}}}(undef, 0)

    for i in axes(p, 1)
        if length(uniq_p) == 0
            push!(uniq_p, p[i])
            push!(uniq_ops, [operators[i]])
            continue
        end

        if !(last(uniq_p) ≈ p[i])
            push!(uniq_p, p[i])
            push!(uniq_ops, [operators[i]])
        else
            push!(uniq_ops[end], operators[i])
        end
    end

    # rescale uniq_p
    for i in eachindex(uniq_p, uniq_ops)
        uniq_p[i] *= length(uniq_ops[i])
    end
    return uniq_ops, uniq_p
end

struct HierarchicalOperatorSampler{N, T, P} <: AbstractOperatorSampler{N, T, P}
    operator_bins::Vector{Vector{NTuple{N, Int}}}
    pvec::P
end

function HierarchicalOperatorSampler(operators::Vector{NTuple{N, Int}}, p::AbstractVector{T}) where {N, T <: Real}
    @assert length(operators) == length(p) "Given vectors must have the same length!"
    operator_bins, p = cluster_probs_vec(operators, p)
    pvec = probability_vector(p)
    return HierarchicalOperatorSampler{N, T, typeof(pvec)}(operator_bins, pvec)
end
length(os::HierarchicalOperatorSampler) = sum(length, os.operator_bins)

function rand(os::HierarchicalOperatorSampler{N})::NTuple{N, Int} where N
    @inbounds ops_list = os.operator_bins[rand(os.pvec)]
    l = rand(1:length(ops_list))

    return @inbounds ops_list[l]
end


text:
function test_belief()
    p = LightDark1D()
    ss0 = MCVI.MCVISubspace([LightDark1DState(0, 4), LightDark1DState(0, 1)], Float64[], Dict{Int64, MCVI.MCVISubspace{LightDark1DState,Int64}}())
    b0 = MCVI.MCVIBelief(ss0, [0.5, 0.5], UInt64(1))
    b1 = MCVI.next(b0, 1, p, MersenneTwister(42))
    @test MCVI.particle(b1.space, 1).y == 5
    b2 = MCVI.next(b0, 5.0, p)
    @test !isnan(b2.weights[1])
    @test b2.weights[1] >= b2.weights[2]
    return true
end


text:
# # Finding balance and variability of constraint-based models

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/@__NAME__.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/@__NAME__.ipynb)

# Here we use [`flux_balance_analysis`](@ref),
# [`flux_variability_analysis`](@ref), and
# [`parsimonious_flux_balance_analysis`](@ref) of `COBREXA.jl` functions to
# analyze a toy model of *E. coli*.

# If it is not already present, download the model.

!isfile("e_coli_core.xml") &&
    download("http://bigg.ucsd.edu/static/models/e_coli_core.xml", "e_coli_core.xml")

using COBREXA

#md # !!! tip "Tip: use `?` to get quick help about functions"
#md #       When you are unsure about how a function works, write `?
#md #       function_name` to see the function reference documentation.

model = load_model("e_coli_core.xml")

# ## Optimization solvers in `COBREXA`
#
# To actually perform any optimization based analysis we need to load an
# optimizer. Any [`JuMP.jl`-supported
# optimizers](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers)
# will work. Here, we will use [`Tulip.jl`](https://github.com/ds4dm/Tulip.jl)
# to optimize linear programs and
# [`OSQP.jl`](https://osqp.org/docs/get_started/julia.html) to optimize quadratic
# programs.

import Pkg
Pkg.add("Tulip")
Pkg.add("OSQP")

using Tulip
using OSQP

# ## Flux balance analysis (FBA)
#
# Most analysis functions come in several variants that produce different types
# of output. All of them usually require a model and `JuMP.jl`-compatible
# optimizer to work in the model.
#
# In the case of FBA, you may choose from these variants (here using the
# `Tulip` optimizer):

vec_soln = flux_balance_analysis_vec(model, Tulip.Optimizer)
#
dict_soln = flux_balance_analysis_dict(model, Tulip.Optimizer)

# ## Modifications

# Often it is desirable to add a slight modififaction to the problem before
# performing analysis, to see e.g. differences of the model behavior caused by
# the change introduced.
#
# `COBREXA.jl` supports several modifications by default, which include
# changing objective sense, optimizer attributes, flux constraints,
# optimization objective, reaction and gene knockouts, and others.

dict_soln = flux_balance_analysis_dict(
    model,
    OSQP.Optimizer;
    modifications = [ # modifications are applied in order
        ## this changes the objective to maximize the biomass production
        change_objective("R_BIOMASS_Ecoli_core_w_GAM"),

        ## this fixes a specific rate of the glucose exchange
        change_constraint("R_EX_glc__D_e", -12, -12),

        ## this knocks out two genes, i.e. constrains their associated reactions to zero.
        knockout(["b0978", "b0734"]), ## the gene IDs are cytochrome oxidase (CYTBD)

        ## ignore the optimizer specified above and change it to Tulip
        change_optimizer(Tulip.Optimizer),

        ## set a custom attribute of the Tulip optimizer (see Tulip docs for more possibilities)
        change_optimizer_attribute("IPM_IterationsLimit", 110),

        ## explicitly tell the optimizer to maximize the new objective
        change_sense(MAX_SENSE),
    ],
)

# ## Flux variability analysis (FVA)

# The default FVA in [`flux_variability_analysis`](@ref) returns maximized and
# minimized reaction fluxes in a matrix. Here we use the dictionary variant in
# flux_variability_analysis_dict, to show how to easily access specific fluxes
# from its results.

fva_mins, fva_maxs = flux_variability_analysis_dict(
    model,
    Tulip.Optimizer;
    bounds = objective_bounds(0.99), # the objective function is allowed to vary by ~1% from the FBA optimum
    modifications = [
        change_optimizer_attribute("IPM_IterationsLimit", 500),
        change_constraint("R_EX_glc__D_e", -10, -10),
        change_constraint("R_EX_o2_e", 0.0, 0.0),
    ],
)

#
fva_maxs["R_EX_ac_e"]["R_EX_ac_e"] # get the maximal acetate exchange flux

# ## Parsimonious flux balance analysis (pFBA)

# Parsimonious flux balance analysis (here in
# [`parsimonious_flux_balance_analysis`](@ref) finds a unique flux solution
# that minimizes the squared sum of fluxes of the system subject, while
# maintaining the same objective value as the flux balance analysis solution.
# Since we are optimizing a quadratic objective, we also need to switch to a
# quadratic optimizer. In this case, OSQP will work. We demonstrate it on the
# dictionary-returning variant of pFBA,
# [`parsimonious_flux_balance_analysis_dict`](@ref):

dict_soln = parsimonious_flux_balance_analysis_dict(
    model,
    OSQP.Optimizer;
    modifications = [
        change_optimizer_attribute("verbose", false), # silence the optimizer (OSQP is very verbose by default)
        change_constraint("R_EX_glc__D_e", -12, -12),
    ],
)

# The function also has the expectable second variant that returns a vector of
# solutions, in [`parsimonious_flux_balance_analysis_vec`](@ref). Here, we
# utilize it to show how to use different optimizers for finding the optimum
# and for solving the quadratic problem. That may be preferable if the
# optimizer qualities differ for the differing tasks. pFBA allows you to
# specify `qp_modifications` that are applied after the original optimum is
# found, and before the quadratic part of the problem solving begins.
vec_soln = parsimonious_flux_balance_analysis_vec(
    model,
    Tulip.Optimizer; # start with Tulip
    modifications = [
        change_constraint("R_EX_glc__D_e", -12, -12),
        change_optimizer_attribute("IPM_IterationsLimit", 500), # we may change Tulip-specific attributes here
    ],
    qp_modifications = [
        change_optimizer(OSQP.Optimizer), # now switch to OSQP (Tulip wouldn't be able to finish the computation)
        change_optimizer_attribute("verbose", false), # and silence it.
    ],
)


text:
using ProbNumDiffEq
using Test
using LinearAlgebra

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems;
importodeproblems();
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_lotkavoltera

h = rand()
σ = rand()

@testset "Test vanilla (ie. non-preconditioned) IBM (d=2,q=2)" begin
    d, q = 2, 2

    A!, Q! = ProbNumDiffEq.vanilla_ibm(d, q)
    Ah = diagm(0 => ones(d * (q + 1)))
    Qh = zeros(d * (q + 1), d * (q + 1))
    A!(Ah, h)
    Q!(Qh, h, σ^2)

    AH_22_IBM = [
        1 h h^2/2 0 0 0
        0 1 h 0 0 0
        0 0 1 0 0 0
        0 0 0 1 h h^2/2
        0 0 0 0 1 h
        0 0 0 0 0 1
    ]
    @test AH_22_IBM ≈ Ah

    QH_22_IBM =
        σ^2 .* [
            h^5/20 h^4/8 h^3/6 0 0 0
            h^4/8 h^3/3 h^2/2 0 0 0
            h^3/6 h^2/2 h 0 0 0
            0 0 0 h^5/20 h^4/8 h^3/6
            0 0 0 h^4/8 h^3/3 h^2/2
            0 0 0 h^3/6 h^2/2 h
        ]
    @test QH_22_IBM ≈ Qh
end

@testset "Test IBM with preconditioning (d=1,q=2)" begin
    d, q = 2, 2

    A, Q = ProbNumDiffEq.ibm(d, q)
    Qh = Q * σ^2

    AH_22_PRE = [
        1 2 1 0 0 0
        0 1 1 0 0 0
        0 0 1 0 0 0
        0 0 0 1 2 1
        0 0 0 0 1 1
        0 0 0 0 0 1
    ]

    QH_22_PRE =
        σ^2 * [
            1/5 1/4 1/3 0 0 0
            1/4 1/3 1/2 0 0 0
            1/3 1/2 1/1 0 0 0
            0 0 0 1/5 1/4 1/3
            0 0 0 1/4 1/3 1/2
            0 0 0 1/3 1/2 1/1
        ]

    @test AH_22_PRE ≈ A
    @test QH_22_PRE ≈ Qh
end

@testset "Verify correct prior dim" begin
    prob = prob_ode_lotkavoltera
    d = length(prob.u0)
    for q in 1:5
        integ = init(prob, EK0(order=q), initialize_derivatives=false)
        @test length(integ.cache.x.μ) == d * (q + 1)
        sol = solve!(integ)
        @test length(integ.cache.x.μ) == d * (q + 1)
        @test length(sol.x_filt[end].μ) == d * (q + 1)
    end
end


text:
# Dynamically keep track of the eigenvalues system of a hermitian, positive
# definite matrix that is growing diagonally.
#   * We need to keep track of its inverse, which is not going to be very stable numerically, but it should do the trick
#   of tracking eigenvalues using the Cauchy Interlace Theorem.
#   * We also want to keep track of it's eigensystem, which is going to serve for the Lanczos Algorithm with selective
#   re-orthogonalizations!
#   * It starts keeping track of the eigensystem when it's been asked to, and it will start doing it.

mutable struct DynamicSymTridiagonal{T<:AbstractFloat}
    alphas::Vector{T}       # Diaognal
    betas::Vector{T}        # Lower & upper Diaognal
    L::Vector{T}            # The lower diagonal of the unit-bidiagonal matrix L
    U::Vector{T}            # Lower diagonal of the upper bi-diagonal matrix U
    k::Int64                # Size of the matrix.

    V::Matrix{T}            # Eigenvectors
    thetas::Vector{T}       # Eigen values.

    function DynamicSymTridiagonal{T}(alpha::T) where {T<:Float64}
        this = new{T}()
        this.alphas = Vector{T}()
        push!(this.alphas, alpha)
        this.betas = Vector{T}()
        this.L = Vector{T}()
        this.U = Vector{T}()
        push!(this.U, alpha)
        this.k = 1
        return this
    end

    function DynamicSymTridiagonal(alpha::AbstractFloat)
        T = typeof(alpha)
        return DynamicSymTridiagonal{T}(alpha)
    end
end

"""
    append an alpha and a beta (diagonal and sub & super diagonal) elemnt to the current sym diagonal matrix.
"""
function (this::DynamicSymTridiagonal{T})(alpha::T, beta::T) where {T<:AbstractFloat}
    push!(this.alphas, alpha)
    push!(this.betas, beta)
    push!(this.U, alpha - beta^2/this.U[end])
    push!(this.L, beta/this.U[end - 1])
    this.k += 1
    return this
end

"""
    Using the Lapack library to get the eigen system of the curent T, and then after
    this, it will keep track of the eigensystem as elements are added to the dynamic matrix.
"""
function InitializeEigenSystem(this::DynamicSymTridiagonal{T}, max_size::Int64) where {T<:AbstractFloat}
    if max_size > 32728 || max_size < 2
        error("Eigen system estimated max size can't be: $(max_size)")
    end
    Trid = GetT(this)
    theta, V = eigen(Trid)
    Vbigger = zeros(T, max_size, max_size)
    Vbigger[1:this.k, 1:this.k] = V
    this.V = Vbigger
    Theta = zeros(T, max_size)
    Theta[1:this.k] = theta
    this.thetas = Theta
    # Sort by absolute values.
    SortedIdx = sortperm(abs.(this.thetas))
    this.thetas = this.thetas[SortedIdx]
    this.V = this.V[:, SortedIdx]
    return this.thetas, this.V
end

"""
    Update the eigensystem of the bigger matrix using the preivous eigen system. Use Power Iterations
    iteratively.
"""
function UpdateEigenSystemPowItr(this::DynamicSymTridiagonal{T}) where {T <: AbstractFloat}
    # if !(isdefined(Base, :this.V))  # cheks if defined first.
    #     error("Must initialize the eigensystem with an estimated size before calling this function. ")
    # end
    k = this.k
    V = this.V
    Θ = this.thetas
    A = this
    if k > size(this.V, 1)
        error("Hasn't implemented this part yet, k exceed the maximal size of eigen system.")
    end
    Avⱼ = zeros(T, k)
    for j in 1: k
        vⱼ = view(V, 1:k, j)
        vⱼ .+= 1e-10*rand(T, k)
        if j >= 2   # Components on to previous ortho eigen vectors removed.
            vⱼ .-= view(V, 1:k, 1:j - 1)*view(V, 1:k, 1:j - 1)'*vⱼ
        end
        for Itr in 1: k^2 + 1e3      # Power iterations
            vⱼᵀvⱼ = dot(vⱼ, vⱼ)
            Apply!(A, vⱼ, Avⱼ)
            if j >= 2  # Components on to previous ortho eigen vectors removed.
                Avⱼ .-= view(V, 1:k, 1:j - 1)*view(V, 1:k, 1:j - 1)'*Avⱼ
            end
            λ̃ = dot(vⱼ, Avⱼ)/vⱼᵀvⱼ  # Rayleigh Quotient.
            ∇r = 2(Avⱼ - λ̃*vⱼ)      # Gradient of Rayleigh Quotient.
            if norm(∇r, Inf) <= 1e-8
                V[1:k, j] .= vⱼ/sqrt(vⱼᵀvⱼ)
                Θ[j] = λ̃
                break               # go to Next Eigenvector.
            end
            vⱼ .= Avⱼ/norm(Avⱼ)
        end
    end
    return
end

"""
    Get the L matrix for this instance.
"""
function GetL(this::DynamicSymTridiagonal)
    return Bidiagonal(fill(1, this.k), convert(Vector{Float64}, this.L), :L)
end

"""
    Get the U matrix for this instance.
"""
function GetU(this::DynamicSymTridiagonal)
    return Bidiagonal(convert(Vector{Float64}, this.U), this.betas, :U)
end

"""
    Get the Tridiagonal Matrix
"""
function GetT(this::DynamicSymTridiagonal)
    return SymTridiagonal(this.alphas, this.betas)
end

"""
    Solve this matrix against another matrix or vector.
"""
function Base.:\(
    this::DynamicSymTridiagonal{T},
    b::Union{Matrix{T}, Vector{T}}
    ) where {T<:AbstractFloat}

    return
end

"""
    Multiply a vector on the left hand size of this dynamically growing matrix
    T.
"""
function Base.:*(this::DynamicSymTridiagonal{T}, b::AbstractArray{T}) where {T <: AbstractFloat}
    @assert length(b) == this.k "Dimension of vector b does match the number of rows in k. Expect $(this.k) but get $(length(b))"
    β = this.betas
    α = this.alphas
    v = similar(b)
    v[1] = α[1]*b[1] + β[1]*b[2]
    for Idx in 2: this.k - 1
        v[Idx] = β[Idx - 1]*b[Idx - 1] + α[Idx]*b[Idx] + β[Idx]*b[Idx + 1]
    end
    v[this.k] = β[this.k - 1]*b[end - 1] + α[this.k]*b[end]
    return v
end

"""
    Memory friendly version of multiplications.
"""
function Apply!(this::DynamicSymTridiagonal{T}, b::AbstractArray{T}, v::Vector{T}) where {T <: AbstractFloat}
    @assert length(b) == this.k "Dimension of vector b does match the number of rows in k. Expect $(this.k) but get $(length(b))"
    @assert length(v) == this.k "Dimension of mutable vector v deosn't match the number of rows in k. Expect $(this.k) but get $(length(v))"
    β = this.betas
    α = this.alphas
    v[1] = α[1]*b[1] + β[1]*b[2]
    for Idx in 2: this.k - 1
        v[Idx] = β[Idx - 1]*b[Idx - 1] + α[Idx]*b[Idx] + β[Idx]*b[Idx + 1]
    end
    v[this.k] = β[this.k - 1]*b[end - 1] + α[this.k]*b[end]
    return v
end

# BASIC TESTING ----------------------------------------------------------------
using LinearAlgebra, Logging, Plots
@info "Basic Testing"
n = 64
A = SymTridiagonal(fill(-2.0, n), fill(1.0,n - 1))
T_dy = DynamicSymTridiagonal(A[1, 1])
InitializeEigenSystem(T_dy, n)
for Idx in 2: n
    T_dy(A[Idx, Idx], A[Idx - 1, Idx])
    @time UpdateEigenSystemPowItr(T_dy)
    # display(T_dy.V)
end
@info "Dynamic T Eigen System"
display(T_dy.V)
L = GetL(T_dy)
U = GetU(T_dy)
b = rand(n)
@assert norm(GetT(T_dy) - L*U) < 1e-10
@assert norm(T_dy*b - GetT(T_dy)*b) < 1e-10


text:
import Base:run

export
run,
Session,
init,
run_profile,
save_profile

function Session(args...;kwargs...)
    config = tf.compat.v1.ConfigProto()
    config.gpu_options.allow_growth=true
    tf.compat.v1.Session(args...;config=config, kwargs...)
end

function Base.:run(o::PyObject, fetches::Union{PyObject, Array{PyObject}, Array{Any}, Tuple}, args::Pair{PyObject, <:Any}...; kwargs...)
    if length(args)>0
        o.run(fetches, feed_dict = Dict(args))
    else
        o.run(fetches; kwargs...)
    end
end

function global_variables_initializer()
    tf.compat.v1.global_variables_initializer()
end

function init(o::PyObject)
    run(o, global_variables_initializer())
end

"""
    run_profile(args...;kwargs...)

Runs the session with tracing information.
"""
function run_profile(args...;kwargs...)
    global run_metadata
    options = tf.compat.v1.RunOptions(trace_level=tf.compat.v1.RunOptions.FULL_TRACE)
    run_metadata = tf.compat.v1.RunMetadata()
    run(args...;options=options, run_metadata=run_metadata, kwargs...)
end

"""
    save_profile(filename::String="default_timeline.json")

Save the timeline information to file `filename`.
- Open Chrome and navigate to chrome://tracing
- Load the timeline file
"""
function save_profile(filename::String="default_timeline.json"; kwargs...)
    timeline = pyimport("tensorflow.python.client.timeline")
    fetched_timeline = timeline.Timeline(run_metadata.step_stats)
    chrome_trace = fetched_timeline.generate_chrome_trace_format(;kwargs...)
    open(filename,"w") do io
        write(io, chrome_trace)
    end
    @info "Timeline information saved in $filename
- Open Chrome and navigate to chrome://tracing
- Load the timeline file"
end


text:
# This script tests the internal functions of eqcond.jl
using Test, OrderedCollections
include(joinpath(dirname(@__FILE__), "../../../src/includeall.jl"))

rp = joinpath(dirname(@__FILE__), "../../reference/models/li2020")

# Set up input arguments and expected behavior
inside_output = load(joinpath(rp, "inside_iteration_loop_output.jld2"))
prepQ̂_input = load(joinpath(rp, "post_processing_input.jld2"))
aug_out = load(joinpath(rp, "augment_variables_out.jld2"), "results")

m = Li2020()
stategrid, funcvar, derivs, endo = initialize!(m)
funcvar[:p] = vec(prepQ̂_input["p_new"])
funcvar[:Q̂] = vec(prepQ̂_input["hat_Q_last"])
funcvar[:xg] = vec(prepQ̂_input["xg_new"])
endo[:κp] = vec(prepQ̂_input["kappap_vec"])
θ = parameters_to_named_tuple(m)
f = eqcond(m)
augment_variables!(m, stategrid, funcvar, derivs, endo)

@testset "Augment variables of Li2020" begin
    @test @test_matrix_approx_eq vec(aug_out["psi_vec"]) endo[:ψ]
    @test @test_matrix_approx_eq vec(aug_out["xK_vec"]) endo[:xK]
    @test @test_matrix_approx_eq vec(aug_out["Lvg"]) endo[:lvg]
    @test @test_matrix_approx_eq vec(aug_out["yK_vec"]) endo[:yK]
    @test @test_matrix_approx_eq vec(aug_out["yg_vec"]) endo[:yg]
    @test @test_matrix_approx_eq vec(aug_out["sigmap_vec"]) endo[:σp]
    @test @test_matrix_approx_eq vec(aug_out["sigmab_vec"]) endo[:σ]
    @test @test_matrix_approx_eq vec(aug_out["sigmah_vec"]) endo[:σh]
    @test @test_matrix_approx_eq vec(aug_out["firesale_jump"]) endo[:firesale_jump]
    @test @test_matrix_approx_eq vec(aug_out["kappad_vec"]) endo[:κd]
    @test @test_matrix_approx_eq vec(aug_out["kappab_vec"]) endo[:κb]
    @test @test_matrix_approx_eq vec(aug_out["kappafs_vec"]) endo[:κfs]
    @test @test_matrix_approx_eq vec(aug_out["kappah_vec"]) endo[:κh]
    @test @test_matrix_approx_eq vec(aug_out["liq_prem_vec"]) endo[:liq_prem]
    @test @test_matrix_approx_eq vec(aug_out["bank_liq_frac"]) endo[:bank_liq_frac]
    @test @test_matrix_approx_eq vec(aug_out["muR_rd_vec"]) endo[:μR_rd]
    @test @test_matrix_approx_eq vec(aug_out["mub_muh_vec"]) endo[:μb_μh]
    @test @test_matrix_approx_eq vec(aug_out["mup_vec"]) endo[:μp]
    @test @test_matrix_approx_eq vec(aug_out["muK_vec"]) endo[:μK]
    @test @test_matrix_approx_eq vec(aug_out["muR_vec"]) endo[:μR]
    @test @test_matrix_approx_eq vec(aug_out["kappaK_vec"]) endo[:κK]
    @test @test_matrix_approx_eq vec(aug_out["rd_vec"]) endo[:rd]
    @test @test_matrix_approx_eq vec(aug_out["rd_vec"]) endo[:rd]
    @test @test_matrix_approx_eq vec(aug_out["muw_vec"]) endo[:μw]
    @test @test_matrix_approx_eq vec(aug_out["sigmaw_vec"]) endo[:σw]
    @test @test_matrix_approx_eq vec(aug_out["kappaw_vec"]) endo[:κw]
    @test @test_matrix_approx_eq vec(aug_out["rg_vec"]) endo[:rg]
    @test @test_matrix_approx_eq vec(aug_out["rh_vec"]) endo[:rh]
    @test @test_matrix_approx_eq vec(aug_out["rf_vec"]) endo[:rf]
end

nothing


text:
export pr_constraint

const PRCache = Dict{Tuple{ProbCircuit, LogicCircuit}, Float64}

# Arthur Choi, Guy Van den Broeck, and Adnan Darwiche. Tractable learning for structured probability
# spaces: A case study in learning preference distributions. In Proceedings of IJCAI, 2015.

"""
Calculate the probability of the logic formula given by LC for the PC
"""
function pr_constraint(pc_node::StructProbCircuit, lc_node, cache::PRCache=PRCache())::Float64

    # TODO require that both circuits have an equal vtree for safety. If they don't, then first convert them to have a vtree
    @assert respects_vtree(lc_node, vtree(pc_node)) "Both circuits do not have an equal vtree"

    # Cache hit
    if (pc_node, lc_node) in keys(cache)
        return cache[pc_node, lc_node]

    # Boundary cases
    elseif isliteralgate(pc_node)
        # Both are literals, just check whether they agrees with each other
        if isliteralgate(lc_node)
            if literal(pc_node) == literal(lc_node)
                return get!(cache, (pc_node, lc_node), 1.0)
            else
                return get!(cache, (pc_node, lc_node), 0.0)
            end
        else
            pr_constraint(pc_node, children(lc_node)[1], cache)
            if length(children(lc_node)) > 1
                pr_constraint(pc_node, children(lc_node)[2], cache)
                return get!(cache, (pc_node, lc_node), 1.0)
            else
                return get!(cache, (pc_node, lc_node),
                    literal(children(lc_node)[1]) == literal(pc_node) ? 1.0 : 0.0)
            end
        end

    # The pc is true
    elseif isliteralgate(children(pc_node)[1])
        theta = exp(pc_node.log_probs[1])
        return get!(cache, (pc_node, lc_node),
            theta * pr_constraint(children(pc_node)[1], lc_node, cache) +
            (1.0 - theta) * pr_constraint(children(pc_node)[2], lc_node, cache))

    # Both pcs are not trivial
    else
        prob = 0.0
        for (prob⋀_node, log_theta) in zip(children(pc_node), pc_node.log_probs)
            p = children(prob⋀_node)[1]
            s = children(prob⋀_node)[2]

            theta = exp(log_theta)
            for lc⋀_node in children(lc_node)
                r = children(lc⋀_node)[1]
                t = children(lc⋀_node)[2]
                prob += theta * pr_constraint(p, r, cache) * pr_constraint(s, t, cache)
            end
        end
        return get!(cache, (pc_node, lc_node), prob)
    end
end


text:
#################### Dependent ####################

const depfxargs = [(:model, Mamba.Model)]

#################### Base Methods ####################

function Base.show(io::IO, d::AbstractDependent)
  msg = string(ifelse(isempty(d.monitor), "An un", "A "),
               "monitored node of type \"", summary(d), "\"\n")
  print(io, msg)
  show(io, d.value)
end

function Base.showall(io::IO, d::AbstractDependent)
  show(io, d)
  print(io, "\nFunction:\n")
  show(io, "text/plain", first(code_typed(d.eval)))
  print(io, "\n\nSource Nodes:\n")
  show(io, d.sources)
  print(io, "\n\nTarget Nodes:\n")
  show(io, d.targets)
end

dims(d::AbstractDependent) = size(d)

function names(d::AbstractDependent)
  names(d, d.symbol)
end

function setmonitor!(d::AbstractDependent, monitor::Bool)
  value = monitor ? Int[0] : Int[]
  setmonitor!(d, value)
end

function setmonitor!(d::AbstractDependent, monitor::Vector{Int})
  values = monitor
  n = length(unlist(d))
  if n > 0 && !isempty(monitor)
    if monitor[1] == 0
      values = collect(1:n)
    elseif minimum(monitor) < 1 || maximum(monitor) > n
      throw(BoundsError())
    end
  end
  d.monitor = values
  d
end

#################### Distribution Fallbacks ####################

unlist(d::AbstractDependent, transform::Bool=false) =
  unlist(d, d.value, transform)

unlist(d::AbstractDependent, x::Real, transform::Bool=false) = [x]

unlist(d::AbstractDependent, x::AbstractArray, transform::Bool=false) = vec(x)

relist(d::AbstractDependent, x::AbstractArray, transform::Bool=false) =
  relistlength(d, x, transform)[1]

logpdf(d::AbstractDependent, transform::Bool=false) = 0.0

logpdf(d::AbstractDependent, x, transform::Bool=false) = 0.0

#################### Logical ####################

@promote_scalarvariate ScalarLogical

#################### Constructors ####################

function Logical(f::Function, monitor::Union{Bool, Vector{Int}}=true)
  value = Float64(NaN)
  fx, src = modelfxsrc(depfxargs, f)
  l = ScalarLogical(value, :nothing, Int[], fx, src, Symbol[])
  setmonitor!(l, monitor)
end

function Logical(d::Integer, f::Function,
                 monitor::Union{Bool, Vector{Int}}=true)
  value = Array{Float64}(fill(0, d)...)
  fx, src = modelfxsrc(depfxargs, f)
  l = ArrayLogical(value, :nothing, Int[], fx, src, Symbol[])
  setmonitor!(l, monitor)
end

#################### Updating ####################

function setinits!(l::AbstractLogical, m::Model, ::Any=nothing)
  l.value = l.eval(m)
  setmonitor!(l, l.monitor)
end

function update!(l::AbstractLogical, m::Model)
  l.value = l.eval(m)
  l
end

#################### Distribution Methods ####################

relistlength(d::ScalarLogical, x::AbstractArray, transform::Bool=false) =
  (x[1], 1)

function relistlength(d::ArrayLogical, x::AbstractArray, transform::Bool=false)
  n = length(d)
  value = reshape(x[1:n], size(d))
  (value, n)
end

#################### Stochastic ####################

#################### Base Methods ####################

@promote_scalarvariate ScalarStochastic

function Base.showall(io::IO, s::AbstractStochastic)
  show(io, s)
  print(io, "\n\nDistribution:\n")
  show(io, s.distr)
  print(io, "\nFunction:\n")
  show(io, "text/plain", first(code_typed(s.eval)))
  print(io, "\n\nSource Nodes:\n")
  show(io, s.sources)
  print(io, "\n\nTarget Nodes:\n")
  show(io, s.targets)
end

#################### Constructors ####################

function Stochastic(f::Function, monitor::Union{Bool, Vector{Int}}=true)
  value = Float64(NaN)
  fx, src = modelfxsrc(depfxargs, f)
  s = ScalarStochastic(value, :nothing, Int[], fx, src, Symbol[],
                       NullUnivariateDistribution())
  setmonitor!(s, monitor)
end

function Stochastic(d::Integer, f::Function,
                    monitor::Union{Bool, Vector{Int}}=true)
  value = Array{Float64}(fill(0, d)...)
  fx, src = modelfxsrc(depfxargs, f)
  s = ArrayStochastic(value, :nothing, Int[], fx, src, Symbol[],
                      NullUnivariateDistribution())
  setmonitor!(s, monitor)
end

#################### Updating ####################

function setinits!(s::ScalarStochastic, m::Model, x::Real)
  s.value = convert(Float64, x)
  s.distr = s.eval(m)
  setmonitor!(s, s.monitor)
end

function setinits!(s::ArrayStochastic, m::Model, x::DenseArray)
  s.value = convert(typeof(s.value), copy(x))
  s.distr = s.eval(m)
  if !isa(s.distr, UnivariateDistribution) && dims(s) != dims(s.distr)
    throw(DimensionMismatch("incompatible distribution for stochastic node"))
  end
  setmonitor!(s, s.monitor)
end

function setinits!(s::AbstractStochastic, m::Model, x)
  throw(ArgumentError("incompatible initial value for node : $(s.symbol)"))
end

function update!(s::AbstractStochastic, m::Model)
  s.distr = s.eval(m)
  s
end

#################### Distribution Methods ####################

function unlist(s::AbstractStochastic, transform::Bool=false)
  unlist(s, s.value, transform)
end

function unlist(s::AbstractStochastic, x::Real, transform::Bool=false)
  unlist(s, [x], transform)
end

function unlist(s::AbstractStochastic, x::AbstractArray, transform::Bool=false)
  transform ? unlist_sub(s.distr, link_sub(s.distr, x)) :
              unlist_sub(s.distr, x)
end

function relist(s::AbstractStochastic, x::AbstractArray, transform::Bool=false)
  relistlength(s, x, transform)[1]
end

function relistlength(s::AbstractStochastic, x::AbstractArray,
                      transform::Bool=false)
  value, n = relistlength_sub(s.distr, s, x)
  (transform ? invlink_sub(s.distr, value) : value, n)
end

function logpdf(s::AbstractStochastic, transform::Bool=false)
  logpdf(s, s.value, transform)
end

function logpdf(s::AbstractStochastic, x, transform::Bool=false)
  logpdf_sub(s.distr, x, transform)
end

rand(s::AbstractStochastic) = rand_sub(s.distr, s.value)


text:
# This file is a part of project JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/NodeNumbering.jl/blob/master/LICENSE

using Documenter
using NodeNumbering

makedocs(
    modules = [NodeNumbering],
    sitename = "NodeNumbering.jl",
    format = :html,
    pages = [
             "Introduction" => "index.md",
             "Theory" => "theory.md",
             "API" => "api.md"
            ])


text:
module MyPackage

export f
export f2

f(x) = 10*x

include("./f2.jl")

end


text:
export MachOSymtabEntry, MachOSymbolRef, MachOSymbols

"""
    MachOSymtabEntry

MachO Symbol Table entry type, internally represented as an `nlist` type.
"""
abstract type MachOSymtabEntry{H <: MachOHandle} <: SymtabEntry{H} end
@io struct nlist{H <: MachOHandle} <: MachOSymtabEntry{H}
    n_strx::UInt32
    n_type::UInt8
    n_sect::UInt8
    n_desc::UInt16
    n_value::UInt32
end

@io struct nlist_64{H <: MachOHandle} <: MachOSymtabEntry{H}
    n_strx::UInt32
    n_type::UInt8
    n_sect::UInt8
    n_desc::UInt16
    n_value::UInt64
end

symbol_name(sym::MachOSymtabEntry) = string("strtab@", sym.n_strx)
symbol_value(sym::MachOSymtabEntry) = sym.n_value
symbol_type(sym::MachOSymtabEntry) = sym.n_type
symbol_section(sym::MachOSymtabEntry) = sym.n_sect
symbol_description(sym::MachOSymtabEntry) = sym.n_desc

function isglobal(sym::MachOSymtabEntry)
    return (symbol_type(sym) & N_EXT) != 0
end
function islocal(sym::MachOSymtabEntry)
    return !isglobal(sym) && (symbol_type(sym) == N_UNDF)
end
function isweak(sym::MachOSymtabEntry)
    return (symbol_description(sym) & (N_WEAK_REF | N_WEAK_DEF)) != 0
end
function isundef(sym::MachOSymtabEntry)
    return (symbol_type(sym) == N_UNDF) ||
           (isglobal(sym) && symbol_section(sym) == NO_SECT)
end

"""
    MachOSymbols

MachO container type for `SymtabEntry` objects.
"""
struct MachOSymbols{H <: MachOHandle} <: Symbols{H}
    cmd::MachOLoadCmdRef{H,MachOSymtabCmd{H}}
end

function MachOSymbols(oh::MachOHandle)
    cmds = find(MachOLoadCmds(oh), [MachOSymtabCmd])
    if isempty(cmds)
        error("Mach-O file does not contain Symtab load commands")
    end
    return MachOSymbols(first(cmds))
end

Symbols(oh::MachOHandle) = MachOSymbols(oh::MachOHandle)
function Symbols(cmd::MachOLoadCmdRef{H,MachOSymtabCmd{H}}) where {H <: MachOHandle}
    return MachOSymbols(cmd)
end

handle(syms::MachOSymbols) = handle(syms.cmd)
endof(syms::MachOSymbols) = symtab_num_symbols(syms.cmd)

symtab_entry_offset(syms::MachOSymbols) = symtab_symbols_offset(syms.cmd)
symtab_entry_size(syms::MachOSymbols) = sizeof(symtab_entry_type(syms))
function symtab_entry_type(syms::MachOSymbols{H}) where {H <: MachOHandle}
    if is64bit(handle(syms))
        return nlist_64{H}
    else
        return nlist{H}
    end
end

function getindex(syms::MachOSymbols, idx)
    # Punt off to `getindex_ref`
    return getindex_ref(
        syms,
        symtab_entry_offset(syms),
        symtab_entry_size(syms),
        symtab_entry_type(syms),
        SymbolRef,
        idx
    )
end

# Add an StrTab() override to be able to load the symbol string table
"""
    StrTab(syms::MachOSymbols)

Given a MachO `Symbols` container, construct a `MachOStrTab` that refers to the
associated symbol string table.
"""
StrTab(syms::MachOSymbols) = StrTab(syms.cmd)

struct MachOSymbolRef{H <: MachOHandle} <: SymbolRef{H}
    syms::MachOSymbols{H}
    entry::MachOSymtabEntry{H}
    idx::UInt32
end

function SymbolRef(syms::MachOSymbols, entry::MachOSymtabEntry, idx)
    return MachOSymbolRef(syms, entry, UInt32(idx))
end

deref(sym::MachOSymbolRef) = sym.entry
Symbols(sym::MachOSymbolRef) = sym.syms
symbol_number(sym::MachOSymbolRef) = sym.idx
function symbol_name(sym::MachOSymbolRef)
    return strtab_lookup(StrTab(Symbols(sym)), sym.entry.n_strx)
end


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

# Schur decomposition
struct Schur{Ty,S<:AbstractMatrix} <: Factorization{Ty}
    T::S
    Z::S
    values::Vector
    Schur{Ty,S}(T::AbstractMatrix{Ty}, Z::AbstractMatrix{Ty}, values::Vector) where {Ty,S} = new(T, Z, values)
end
Schur(T::AbstractMatrix{Ty}, Z::AbstractMatrix{Ty}, values::Vector) where {Ty} = Schur{Ty, typeof(T)}(T, Z, values)

"""
    schurfact!(A::StridedMatrix) -> F::Schur

Same as [`schurfact`](@ref) but uses the input argument as workspace.
"""
schurfact!(A::StridedMatrix{<:BlasFloat}) = Schur(LinAlg.LAPACK.gees!('V', A)...)

"""
    schurfact(A::StridedMatrix) -> F::Schur

Computes the Schur factorization of the matrix `A`. The (quasi) triangular Schur factor can
be obtained from the `Schur` object `F` with either `F[:Schur]` or `F[:T]` and the
orthogonal/unitary Schur vectors can be obtained with `F[:vectors]` or `F[:Z]` such that
`A = F[:vectors]*F[:Schur]*F[:vectors]'`. The eigenvalues of `A` can be obtained with `F[:values]`.

# Examples
```jldoctest
julia> A = [-2. 1. 3.; 2. 1. -1.; -7. 2. 7.]
3×3 Array{Float64,2}:
 -2.0  1.0   3.0
  2.0  1.0  -1.0
 -7.0  2.0   7.0

julia> F = schurfact(A)
Base.LinAlg.Schur{Float64,Array{Float64,2}} with factors T and Z:
[2.0 0.801792 6.63509; -8.55988e-11 2.0 8.08286; 0.0 0.0 1.99999]
[0.577351 0.154299 -0.801784; 0.577346 -0.77152 0.267262; 0.577354 0.617211 0.534522]
and values:
Complex{Float64}[2.0+8.28447e-6im, 2.0-8.28447e-6im, 1.99999+0.0im]

julia> F[:vectors] * F[:Schur] * F[:vectors]'
3×3 Array{Float64,2}:
 -2.0  1.0   3.0
  2.0  1.0  -1.0
 -7.0  2.0   7.0
```
"""
schurfact(A::StridedMatrix{<:BlasFloat}) = schurfact!(copy(A))
function schurfact(A::StridedMatrix{T}) where T
    S = promote_type(Float32, typeof(one(T)/norm(one(T))))
    return schurfact!(copy_oftype(A, S))
end

function getindex(F::Schur, d::Symbol)
    if d == :T || d == :Schur
        return F.T
    elseif d == :Z || d == :vectors
        return F.Z
    elseif d == :values
        return F.values
    else
        throw(KeyError(d))
    end
end

function show(io::IO, F::Schur)
    println(io, "$(typeof(F)) with factors T and Z:")
    show(io, F[:T])
    println(io)
    show(io, F[:Z])
    println(io)
    println(io, "and values:")
    show(io, F[:values])
end

"""
    schur(A::StridedMatrix) -> T::Matrix, Z::Matrix, λ::Vector

Computes the Schur factorization of the matrix `A`. The methods return the (quasi)
triangular Schur factor `T` and the orthogonal/unitary Schur vectors `Z` such that
`A = Z*T*Z'`. The eigenvalues of `A` are returned in the vector `λ`.

See [`schurfact`](@ref).

# Examples
```jldoctest
julia> A = [-2. 1. 3.; 2. 1. -1.; -7. 2. 7.]
3×3 Array{Float64,2}:
 -2.0  1.0   3.0
  2.0  1.0  -1.0
 -7.0  2.0   7.0

julia> T, Z, lambda = schur(A)
([2.0 0.801792 6.63509; -8.55988e-11 2.0 8.08286; 0.0 0.0 1.99999], [0.577351 0.154299 -0.801784; 0.577346 -0.77152 0.267262; 0.577354 0.617211 0.534522], Complex{Float64}[2.0+8.28447e-6im, 2.0-8.28447e-6im, 1.99999+0.0im])

julia> Z * T * Z'
3×3 Array{Float64,2}:
 -2.0  1.0   3.0
  2.0  1.0  -1.0
 -7.0  2.0   7.0
```
"""
function schur(A::StridedMatrix)
    SchurF = schurfact(A)
    SchurF[:T], SchurF[:Z], SchurF[:values]
end
schur(A::Symmetric) = schur(full(A))
schur(A::Hermitian) = schur(full(A))
schur(A::UpperTriangular) = schur(full(A))
schur(A::LowerTriangular) = schur(full(A))
schur(A::Tridiagonal) = schur(Matrix(A))

"""
    ordschur!(F::Schur, select::Union{Vector{Bool},BitVector}) -> F::Schur

Same as [`ordschur`](@ref) but overwrites the factorization `F`.
"""
function ordschur!(schur::Schur, select::Union{Vector{Bool},BitVector})
    _, _, vals = ordschur!(schur.T, schur.Z, select)
    schur[:values][:] = vals
    return schur
end

"""
    ordschur(F::Schur, select::Union{Vector{Bool},BitVector}) -> F::Schur

Reorders the Schur factorization `F` of a matrix `A = Z*T*Z'` according to the logical array
`select` returning the reordered factorization `F` object. The selected eigenvalues appear
in the leading diagonal of `F[:Schur]` and the corresponding leading columns of
`F[:vectors]` form an orthogonal/unitary basis of the corresponding right invariant
subspace. In the real case, a complex conjugate pair of eigenvalues must be either both
included or both excluded via `select`.
"""
ordschur(schur::Schur, select::Union{Vector{Bool},BitVector}) =
    Schur(ordschur(schur.T, schur.Z, select)...)

"""
    ordschur!(T::StridedMatrix, Z::StridedMatrix, select::Union{Vector{Bool},BitVector}) -> T::StridedMatrix, Z::StridedMatrix, λ::Vector

Same as [`ordschur`](@ref) but overwrites the input arguments.
"""
ordschur!(T::StridedMatrix{Ty}, Z::StridedMatrix{Ty}, select::Union{Vector{Bool},BitVector}) where {Ty<:BlasFloat} =
    LinAlg.LAPACK.trsen!(convert(Vector{BlasInt}, select), T, Z)[1:3]

"""
    ordschur(T::StridedMatrix, Z::StridedMatrix, select::Union{Vector{Bool},BitVector}) -> T::StridedMatrix, Z::StridedMatrix, λ::Vector

Reorders the Schur factorization of a real matrix `A = Z*T*Z'` according to the logical
array `select` returning the reordered matrices `T` and `Z` as well as the vector of
eigenvalues `λ`. The selected eigenvalues appear in the leading diagonal of `T` and the
corresponding leading columns of `Z` form an orthogonal/unitary basis of the corresponding
right invariant subspace. In the real case, a complex conjugate pair of eigenvalues must be
either both included or both excluded via `select`.
"""
ordschur(T::StridedMatrix{Ty}, Z::StridedMatrix{Ty}, select::Union{Vector{Bool},BitVector}) where {Ty<:BlasFloat} =
    ordschur!(copy(T), copy(Z), select)

struct GeneralizedSchur{Ty,M<:AbstractMatrix} <: Factorization{Ty}
    S::M
    T::M
    alpha::Vector
    beta::Vector{Ty}
    Q::M
    Z::M
    function GeneralizedSchur{Ty,M}(S::AbstractMatrix{Ty}, T::AbstractMatrix{Ty}, alpha::Vector,
                                    beta::Vector{Ty}, Q::AbstractMatrix{Ty}, Z::AbstractMatrix{Ty}) where {Ty,M}
        new(S, T, alpha, beta, Q, Z)
    end
end
function GeneralizedSchur(S::AbstractMatrix{Ty}, T::AbstractMatrix{Ty}, alpha::Vector,
                          beta::Vector{Ty}, Q::AbstractMatrix{Ty}, Z::AbstractMatrix{Ty}) where Ty
    GeneralizedSchur{Ty, typeof(S)}(S, T, alpha, beta, Q, Z)
end

"""
    schurfact!(A::StridedMatrix, B::StridedMatrix) -> F::GeneralizedSchur

Same as [`schurfact`](@ref) but uses the input matrices `A` and `B` as workspace.
"""
schurfact!(A::StridedMatrix{T}, B::StridedMatrix{T}) where {T<:BlasFloat} =
    GeneralizedSchur(LinAlg.LAPACK.gges!('V', 'V', A, B)...)

"""
    schurfact(A::StridedMatrix, B::StridedMatrix) -> F::GeneralizedSchur

Computes the Generalized Schur (or QZ) factorization of the matrices `A` and `B`. The
(quasi) triangular Schur factors can be obtained from the `Schur` object `F` with `F[:S]`
and `F[:T]`, the left unitary/orthogonal Schur vectors can be obtained with `F[:left]` or
`F[:Q]` and the right unitary/orthogonal Schur vectors can be obtained with `F[:right]` or
`F[:Z]` such that `A=F[:left]*F[:S]*F[:right]'` and `B=F[:left]*F[:T]*F[:right]'`. The
generalized eigenvalues of `A` and `B` can be obtained with `F[:alpha]./F[:beta]`.
"""
schurfact(A::StridedMatrix{T},B::StridedMatrix{T}) where {T<:BlasFloat} = schurfact!(copy(A),copy(B))
function schurfact(A::StridedMatrix{TA}, B::StridedMatrix{TB}) where {TA,TB}
    S = promote_type(Float32, typeof(one(TA)/norm(one(TA))), TB)
    return schurfact!(copy_oftype(A, S), copy_oftype(B, S))
end

"""
    ordschur!(F::GeneralizedSchur, select::Union{Vector{Bool},BitVector}) -> F::GeneralizedSchur

Same as `ordschur` but overwrites the factorization `F`.
"""
function ordschur!(gschur::GeneralizedSchur, select::Union{Vector{Bool},BitVector})
    _, _, α, β, _, _ = ordschur!(gschur.S, gschur.T, gschur.Q, gschur.Z, select)
    gschur[:alpha][:] = α
    gschur[:beta][:] = β
    return gschur
end

"""
    ordschur(F::GeneralizedSchur, select::Union{Vector{Bool},BitVector}) -> F::GeneralizedSchur

Reorders the Generalized Schur factorization `F` of a matrix pair `(A, B) = (Q*S*Z', Q*T*Z')`
according to the logical array `select` and returns a GeneralizedSchur object `F`. The
selected eigenvalues appear in the leading diagonal of both `F[:S]` and `F[:T]`, and the
left and right orthogonal/unitary Schur vectors are also reordered such that
`(A, B) = F[:Q]*(F[:S], F[:T])*F[:Z]'` still holds and the generalized eigenvalues of `A`
and `B` can still be obtained with `F[:alpha]./F[:beta]`.
"""
ordschur(gschur::GeneralizedSchur, select::Union{Vector{Bool},BitVector}) =
    GeneralizedSchur(ordschur(gschur.S, gschur.T, gschur.Q, gschur.Z, select)...)

"""
    ordschur!(S::StridedMatrix, T::StridedMatrix, Q::StridedMatrix, Z::StridedMatrix, select) -> S::StridedMatrix, T::StridedMatrix, Q::StridedMatrix, Z::StridedMatrix, α::Vector, β::Vector

Same as [`ordschur`](@ref) but overwrites the factorization the input arguments.
"""
ordschur!(S::StridedMatrix{Ty}, T::StridedMatrix{Ty}, Q::StridedMatrix{Ty},
    Z::StridedMatrix{Ty}, select::Union{Vector{Bool},BitVector}) where {Ty<:BlasFloat} =
        LinAlg.LAPACK.tgsen!(convert(Vector{BlasInt}, select), S, T, Q, Z)

"""
    ordschur(S::StridedMatrix, T::StridedMatrix, Q::StridedMatrix, Z::StridedMatrix, select) -> S::StridedMatrix, T::StridedMatrix, Q::StridedMatrix, Z::StridedMatrix, α::Vector, β::Vector

Reorders the Generalized Schur factorization of a matrix pair `(A, B) = (Q*S*Z', Q*T*Z')`
according to the logical array `select` and returns the matrices `S`, `T`, `Q`, `Z` and
vectors `α` and `β`.  The selected eigenvalues appear in the leading diagonal of both `S`
and `T`, and the left and right unitary/orthogonal Schur vectors are also reordered such
that `(A, B) = Q*(S, T)*Z'` still holds and the generalized eigenvalues of `A` and `B` can
still be obtained with `α./β`.
"""
ordschur(S::StridedMatrix{Ty}, T::StridedMatrix{Ty}, Q::StridedMatrix{Ty},
    Z::StridedMatrix{Ty}, select::Union{Vector{Bool},BitVector}) where {Ty<:BlasFloat} =
        ordschur!(copy(S), copy(T), copy(Q), copy(Z), select)

function getindex(F::GeneralizedSchur, d::Symbol)
    if d == :S
        return F.S
    elseif d == :T
        return F.T
    elseif d == :alpha
        return F.alpha
    elseif d == :beta
        return F.beta
    elseif d == :values
        return F.alpha./F.beta
    elseif d == :Q || d == :left
        return F.Q
    elseif d == :Z || d == :right
        return F.Z
    else
        throw(KeyError(d))
    end
end

"""
    schur(A::StridedMatrix, B::StridedMatrix) -> S::StridedMatrix, T::StridedMatrix, Q::StridedMatrix, Z::StridedMatrix, α::Vector, β::Vector

See [`schurfact`](@ref).
"""
function schur(A::StridedMatrix, B::StridedMatrix)
    SchurF = schurfact(A, B)
    SchurF[:S], SchurF[:T], SchurF[:Q], SchurF[:Z], SchurF[:alpha], SchurF[:beta]
end

# Conversion
convert(::Type{AbstractMatrix}, F::Schur) = (F.Z * F.T) * F.Z'
convert(::Type{AbstractArray}, F::Schur) = convert(AbstractMatrix, F)
convert(::Type{Matrix}, F::Schur) = convert(Array, convert(AbstractArray, F))
convert(::Type{Array}, F::Schur) = convert(Matrix, F)
full(F::Schur) = convert(AbstractArray, F)

copy(F::Schur) = Schur(copy(F.T), copy(F.Z), copy(F.values))
copy(F::GeneralizedSchur) = GeneralizedSchur(copy(F.S), copy(F.T), copy(F.alpha), copy(F.beta), copy(F.Q), copy(F.Z))


text:
using CSV, DataFrames, Printf, Random

include("gen_confus.jl")

function gen_dataset(seed = 123, folder="dataset")
  info_filenames = String[]
  info_rows = Int[]
  info_elems = Int[]
  count = 0
  if isdir(folder)
    error("dataset folder exists. Remove it")
  else
    mkdir(folder)
  end
  Random.seed!(seed)
  for col_max = round.(Int, 10 .^ (0.5:0.5:4.5))
    for n_rows = 50:50:500
      for dense in [true, false]
        cm = dense ? col_max : 10col_max
        count += 1
        filename = @sprintf("dataset/file-%04d.txt", count)
        n_elem = gen_confus(filename=filename, n_rows=n_rows, col_max=cm, dense=dense)
        println("$count rows=$n_rows n_elem=$n_elem cm=$cm dense=$dense")
        push!(info_filenames, filename)
        push!(info_rows, n_rows)
        push!(info_elems, n_elem)
      end
    end
  end

  DataFrame(:filename => info_filenames, :nrows => info_rows, :nelements => info_elems) |> CSV.write("dataset/info.csv")
end


text:
using Requests
import Requests: get, post, put, delete, options
using CSV
using DataFrames
using WriteVTK

query = Dict{String, String}(
	"batch" => "l",
	"COMMAND" => "499",
	"CENTER" => "500@10",
	"MAKE_EPHEM" => "YES",
	"TABLE_TYPE" => "VECTORS",
	"START_TIME" => "2017-02-15",
	"STOP_TIME" => "2018-02-15",
	"STEP_SIZE" => "1d",
	"OUT_UNITS" => "KM-S",
	"REF_PLANE" => "ECLIPTIC",
	"REF_SYSTEM" => "J2000",
	"VECT_CORR" => "NONE",
	"VEC_LABELS" => "YES",
	"VEC_DELTA_T" => "NO",
	"CSV_FORMAT" => "YES",
	"OBJ_DATA" => "NO",
	"VEC_TABLE" => "3")

function horizon_data(bodies::Dict{String, String}; url="http://ssd.jpl.nasa.gov/horizons_batch.cgi")
	for (i, k) in enumerate(keys(bodies))
		#download data from horizons
		query["COMMAND"] = k
		response = get(url, query = query)
		println(response)
		str = readstring(response)

		#filter out data
		m = match(r"\$\$SOE\n(?<csv_data>[\s-\S]*)\n\$\$EOE", str)
		data = m[:csv_data]
		data = replace(data, ",\n", "\n")
		df = CSV.read(IOBuffer(data); types=Dict(2=>String), header=["date1", "date2", "X", "Y", "Z", "VX", "VY", "VZ", "LT", "RG", "RR"], nullable=false)
		CSV.write("$(bodies[k]).csv", df)
	end
end

function csv_to_vtu(bodies::Dict{String,String}; file_name="my_vtu_file")
	vectors = Array{DataFrame}(bodies.count)

	for (i, k) in enumerate(keys(bodies))
		vectors[i] = CSV.read("$(bodies[k]).csv", nullable=false)
	end

	cells = WriteVTK.MeshCell[MeshCell(VTKCellTypes.VTK_VERTEX, [i]) for i = 1:length(vectors)]
	for i = 1:size(vectors[1])[1]
		points = Array{Float64}(3, length(vectors))
		for j = 1:length(vectors)
			points[1,j] = vectors[j][:X][i]
			points[2,j] = vectors[j][:Y][i]
			points[3,j] = vectors[j][:Z][i]
		end
		outfile = vtk_grid(string(file_name,i), points, cells) do vtk
			vtk_cell_data(vtk, [1 for i = 1:length(vectors)], "Mass")
		end
	end
end

major_bodies = Dict{String,String}(
	"10" => "Sun",
	"199" => "Mercury",
	"299" => "Venus",
	#"301" => "Moon",
	"399" => "Earth",
	#"401" => "Phobos",
	#"402" => "Deimos",
	"499" => "Mars",
	"599" => "Jupiter",
	"699" => "Saturn",
	"799" => "Uranus",
	"800" => "Neptune")

horizon_data(major_bodies)
csv_to_vtu(major_bodies, file_name="major_bodies")