text:
#DIRECTED NETWORKS------------
# NOT UPDATED TO LAST VERSION, NOT WORKING!
# --------------------------Sender Receiver Effect

struct  SdErgmDirBin1 <: SdErgm
     """ A Logistic Gas model for Directed binary networks and probability depending only
             on time varying parameters.
             ''p^{(t)}_{ij}  = logistic( theta^{(t)}_i + theta^{(t)}_j)''
             Bin1 Model with one unconditional mean per node and one gas parameter per group of nodes
             # Model Definition and parameters conversion functions
         """
     # node - specific W and group Specific (A B)
     obsT:: Array{<:Real,2} # In degrees and out Degrees for each t . The first dimension is always time. Tx2N
     Par::Array{Array{<:Real,1},1} #[[InW;outW], Agroups , Bgroups]  Group specific parameters
     #Nodes in a group have the same GAS parameters, but different fitnesses

     # each node belongs to a W group, each W group is assigned 2 AB groups, in order
     #to allow different parameters for in and our. The AB groups In and Out can be the same

     groupsInds::Array{Array{Int,1},1} #[[Par_2_groupsW, Wgroups_2_groupsAB ]  (_2_ stands for "assignement to" )
     # In and out W parameters are assumed to be  different, while In and Out AB
     # are be equal. Moreover, I assume that nodes with the same Wgroup have the same ABgroups
     scoreScalingType::String # String that specifies the rescaling of the score. For a list of possible
     # choices see function scalingMatGas
  end

#inizializza senza specificare assenza di scaling
SdErgmDirBin1( obsT , Par  , groupsInds) =
        SdErgmDirBin1( obsT,Par, groupsInds,"")
#Initialize by observations only
SdErgmDirBin1(obsT:: Array{<:Real,2}) =   (N = round(length(obsT[:,1])/2);
                                                SdErgmDirBin1(obsT,
                                                    [ zeros(Int(2N)) ,0.9  * ones(1), 0.01 * ones(1)],
                                                    [Int.(1:2N),Int.(ones(Int(2N)))],"") )

SdErgmDirBin1(obsT:: Array{<:Real,2}, scoreScalingType::String, group_string::String) =
    (N = round(length(obsT[:,1])/2);
    if group_string == "ONE_PAR_EACH"
        return SdErgmDirBin1(obsT,
            [ zeros(Int(2N)) ,0.9  * ones(1), 0.01 * ones(1)],
            [Int.(1:2N), Int.(1:2N)], scoreScalingType)
    elseif group_string == "ONE_PAR_ALL"
        return SdErgmDirBin1(obsT,
            [ zeros(Int(2N)) ,0.9  * ones(1), 0.01 * ones(1)],
            [Int.(1:2N),Int.(ones(Int(2N)))],scoreScalingType)
    end)
SdErgmDirBin1(obsT:: Array{<:Real,2},scoreScalingType::String) =   (N = round(length(obsT[:,1])/2);
                                                SdErgmDirBin1(obsT,
                                                    [ zeros(Int(2N)) ,0.9  * ones(1), 0.01 * ones(1)],
                                                    [Int.(1:2N),Int.(ones(Int(2N)))],scoreScalingType) )
fooSdErgmDirBin1 = SdErgmDirBin1(ones(100,20))

# Relations between Static and Dynamical models: conventions on storage for
# parameters and observations
StaModType(Model::SdErgmDirBin1) = StaticNets.fooErgmDirBin1# to be substituted with a conversion mechanism
linSpacedPar(Model::SdErgmDirBin1,Nnodes::Int;NgroupsW = Nnodes, deltaN::Int=3,graphConstr = true) =
        StaticNets.linSpacedPar(StaModType(Model),Nnodes;Ngroups = NgroupsW,deltaN=deltaN,graphConstr =graphConstr);

# options and conversions of parameters for optimization
function setOptionsOptim(Model::SdErgmDirBin1)
    "Set the options for the optimization required in the estimation of the model.
    For the optimization use the Optim package."
    tol = eps()*1000
    opt = Optim.Options(  g_tol = tol,
                     x_tol = tol,
                     f_tol = tol,
                     iterations = 500,
                     show_trace = true,#false,#
                     show_every=1)#,
    #                 extended_trace = true)

    algo = Newton(; alphaguess = LineSearches.InitialHagerZhang(),
                     linesearch = LineSearches.HagerZhang())
    algo = NewtonTrustRegion(; initial_delta = 0.1,
                    delta_hat = 0.8,
                    eta = 0.1,
                    rho_lower = 0.25,
                    rho_upper = 0.75)# Newton(;alphaguess = Optim.LineSearches.InitialStatic(),
       return opt, algo
 end

function array_2_vec_all_par(Model::SdErgmDirBin1,ArrayGasPar::Array{<:Array{Float64,1},1})
         VecGasPar = [ArrayGasPar[1];ArrayGasPar[2];ArrayGasPar[3]]
         return VecGasPar
     end

function NumberOfGroups(Model::SdErgmDirBin1,groupsInds::Array{<:Array{<:Real,1},1})
    NGW = length(unique(groupsInds[1]))
    GBA = length( unique( groupsInds[2][.!(groupsInds[2].==0)] ) )
     return NGW, GBA
 end
function NumberOfGroupsAndABindNodes(Model::SdErgmDirBin1,groupsInds::Array{<:Array{<:Real,1},1})
    NGW = length(unique(groupsInds[1]))
    GBA = length( unique( [groupsInds[2][.!(groupsInds[2].==0)]  ]) )
    ABgroupsIndNodesIO = groupsInds[2][groupsInds[1]]
    indTvNodesIO =   (.!(ABgroupsIndNodesIO .==0))#if at least one group ind is zero then both in and out tv par are constant
    return NGW, GBA, ABgroupsIndNodesIO, indTvNodesIO
 end
 function NumberOfGroupsAndABindNodes(Model::SdErgmDirBin1,groupsInds::Array{Array{Int,1},1})
     NGW = length(unique(groupsInds[1]))
     GBA = length( unique( groupsInds[2][.!(groupsInds[2].==0)] ) )
     ABgroupsIndNodesIO = groupsInds[2][groupsInds[1]]
     indTvNodesIO =   (.!(ABgroupsIndNodesIO .==0))#if at least one group ind is zero then both in and out tv par are constant
     return NGW, GBA, ABgroupsIndNodesIO, indTvNodesIO
  end
function vec_2_array_all_par(Model::SdErgmDirBin1,VecGasPar::Array{<:Real,1};groupsInds::Array{<:Array{<:Real,1},1} = Model.groupsInds)
     println(size(VecGasPar))
     NGW, GBA = NumberOfGroups(Model,groupsInds)
     ArrayGasPar =  [VecGasPar[1:NGW],ones(GBA).*VecGasPar[NGW+1:NGW+GBA],ones(GBA).*VecGasPar[NGW+GBA+1:NGW + 2GBA]]
     return ArrayGasPar
     end
function restrictGasPar(Model::SdErgmDirBin1,vecUnGasPar::Array{<:Real,1})
     "From the Unrestricted values of the parameters return the restricted ones
     takes as inputs a vector of parameters that can take any value in R and returns
     the appropriately contrained values, e.g. the coefficient of autoregressive
     component in the GAS update rule has to be positive and <1."
     groupsInds = Model.groupsInds
     NGW, GBA= NumberOfGroups(Model,groupsInds)
     W_un = vecUnGasPar[1:2NGW]
     diag_B_Un = vecUnGasPar[2NGW+1:2NGW+GBA]
     diag_A_Un = vecUnGasPar[2NGW+ GBA+1:2NGW+2GBA]
     W_Re = W_un
     diag_A_Re = exp.(diag_A_Un)
     diag_B_Re = 1 ./ (1+exp.(.-diag_B_Un))
     vecRePar =  [W_Re; diag_B_Re; diag_A_Re]
     return vecRePar
     end
function unRestrictGasPar( Model::SdErgmDirBin1,vecReGasPar::Array{<:Real,1})
     "From the restricted values of the parameters return the unrestricted ones
     takes as inputs a vector of parameters that can take any value in R and returns
     the appropriately contrained values, i.e. the coefficient of autoregressive
     component in the GAS update rule has to be positive and <1."
     groupsInds = Model.groupsInds
     NGW, GBA= NumberOfGroups(Model,groupsInds)
     W_Re = vecReGasPar[1:2NGW]
     diag_B_Re = vecReGasPar[2NGW+1:2NGW+GBA]
     diag_A_Re = vecReGasPar[2NGW+GBA+1:end]

     W_Un = W_Re
     diag_B_Un = log.(diag_B_Re ./ (1 .- diag_B_Re ))
     diag_A_Un = log.(diag_A_Re)

     vecUnPar =  [W_Re;diag_B_Un; diag_A_Un]
     return vecUnPar
     end

#Gas Filter Functions
identify(Model::SdErgmDirBin1,parIO::Array{Tp,1} where Tp<:Real; idType = "equalIOsums")=
 StaticNets.identify(StaticNets.fooErgmDirBin1,parIO;  idType =idType )

function scalingMatGas(Model::SdErgmDirBin1,expMat::Array{<:Real,2};
                        I_tm1:: Union{UniformScaling, Matrix}=UniformScaling(2))
    "Return the matrix required for the scaling of the score, given the expected
     matrix and the Scaling matrix at previous time. "
    if Model.scoreScalingType == ""
        scalingMat = 1#
    elseif Model.scoreScalingType == "FISHER-EWMA"

        # λ = 0.8
        #
        # I = expMat.*(1-expMat)
        # diagI = sum(I,2)
        # [I[i,i] = diagI[i] for i=1:length(diagI) ]
        # I_t =  λ*I + (1-λ) *I_tm1
        # scalingMat = sqrt(I) ##
    elseif  Model.scoreScalingType == "FISHER-DIAG"
        # display(expMat)
         I = expMat.*(1 .- expMat)
        # scalingMat = zeros(Real,2 .* size(expMat))
        diagScalIn = sqrt.(sumSq(I, 2))
        N = length(diagScalIn)
        # [scalingMat[i,i] = diagScalIn[i] for i=1:N ]
        diagScalOut = sqrt.(sumSq(I, 1))
        # display(diagScalIn)
        # display(diagScalOut)

        # [scalingMat[N+i,N+i] = diagScalOut[i] for i=1:N ]
        scalingMat = Diagonal(vcat(diagScalIn, diagScalOut))
        #print(size(scalingMat))
    end
    return scalingMat
 end
function predict_score_driven_par( Model::SdErgmDirBin1,N::Int,degsIO_t::Array{<:Real,1},
                         ftotIO_t::Array{<:Real,1},I_tm1:: Union{UniformScaling, Matrix},
                         indTvNodesIO::Union{BitArray{1},Array{Bool}},
                         WgasIO::Array{<:Real,1},BgasIO::Array{<:Real,1},AgasIO::Array{<:Real,1})
     #= likelihood and gradients depend on all the parameters (ftot_t), but
     only the time vaying ones (f_t) are to be updated=#

     fIO_t = ftotIO_t[indTvNodesIO] #Time varying fitnesses
     thetas_mat_t_exp, exp_mat_t = StaticNets.expMatrix2(StaticNets.fooErgmDirBin1,ftotIO_t)
     exp_degIO_t = [sumSq(exp_mat_t,2);sumSq(exp_mat_t,1)]
         #The following part is a faster version of logLikelihood_t()
     loglike_t = sum(ftotIO_t .* degsIO_t) -  sum(log.(1 .+ thetas_mat_t_exp))
     # compute the score of observations at time t wrt parameters at tm1
     gradIO_t = degsIO_t  .-  exp_degIO_t # only for TV paramaters
     # updating direction without rescaling (for the moment)
     I_t = scalingMatGas(Model,exp_mat_t;I_tm1=I_tm1)
     # display(diag(I_t))
     # display(size(gradIO_t))
     # display(size(exp_mat_t ))
      sIO_t =   I_t\gradIO_t
      sIO_t[gradIO_t.==0] .= 0

     # GAS update for next iteration
     fIO_tp1 = WgasIO .+ BgasIO .* fIO_t .+ AgasIO .* sIO_t[indTvNodesIO]
     fIO_t = fIO_tp1
     ftotIO_tp1 = copy(ftotIO_t)
     ftotIO_tp1 = fIO_t #of all parameters udate dynamic ones with GAS
     ftotIO_tp1 =identify(Model,ftotIO_tp1 )
     return ftotIO_tp1,loglike_t,gradIO_t
  end
function score_driven_filter_or_dgp( Model::SdErgmDirBin1,
                                vResGasPar::Array{<:Real,1};vConstPar::Array{<:Real,1} = zeros(Real,2),
                                obsT::Array{<:Real,2}=Model.obsT, ftotIO_0::Array{<:Real,1} = zeros(2),
                                groupsInds::Array{<:Array{<:Real,1},1}=Model.groupsInds,
                                dgpNT = (0,0))
    """GAS Filter the Dynamic Fitnesses from the Observed degrees, given the GAS parameters
     given T observations for the degrees in TxN vector degsT
     """
     #se parametri per dgp sono dati allora campiona il dgp, altrimenti filtra
    sum(dgpNT) == 0   ?     dgp = false : dgp = true
    if dgp
        N,T = dgpNT
        Y_T = zeros(Int8,N,N,T)
    else
        N2,T = size(obsT);N = round(Int,N2/2)
    end

    NGW,GBA,ABgroupsIndNodesIO,indTvNodesIO = NumberOfGroupsAndABindNodes(Model,groupsInds)

    # ABgroups Inds for each par (In and Out). IO means that we have a 2N vector
    # of indices: one index for each In par and below one index for each Out par
    #println(NumberOfGroupsAndABindNodes(Model,groupsInds))

    # Organize parameters of the GAS update equation
    WGroupsIO = vResGasPar[1:NGW]
    #    StaticNets.expMatrix2(StaticNets.fooErgmDirBin1,WGroupsIO )
    W_allIO = WGroupsIO[groupsInds[1]]

    BgasGroups  = vResGasPar[NGW+1:NGW+GBA]
    AgasGroups  = vResGasPar[NGW+GBA+1:NGW+2GBA]
    #distribute nodes among groups with equal parameters
    AgasIO = AgasGroups[ABgroupsIndNodesIO[indTvNodesIO]]
    BgasIO = BgasGroups[ABgroupsIndNodesIO[indTvNodesIO]]

    WgasIO   = W_allIO[indTvNodesIO]

    if all(ftotIO_0.==0)
        # start values equal the unconditional mean,
        # but  constant ones remain equal to the unconditional mean, hence initialize as
        ftotIO_0  =  WgasIO ./ (1 .- BgasIO)
    end

    ftotIO_t =  identify(Model, ftotIO_0)
    I_tm1 = UniformScaling(N)
    loglike = 0

    fIOVecT = ftotIO_t
    for t=1:T-1

        if dgp
            expMat_t = StaticNets.expMatrix(StaModType(Model),ftotIO_t )# exp.(ftot_t)
            Y_T[:,:,t] = StaticNets.samplSingMatCan(StaModType(Model),expMat_t)
            degsIO_t = [dropdims(sum(Y_T[:,:,t],dims = 2),dims = 2); dropdims(sum(Y_T[:,:,t],dims = 1),dims = 1)]
        else
            degsIO_t = obsT[:,t] # vector of in and out degrees
        end
        #t==150  ?     println(WgasIO) : ()
        ftotIO_tp1,loglike_t = predict_score_driven_par(Model,N,degsIO_t,ftotIO_t,I_tm1,
                                            indTvNodesIO,WgasIO,BgasIO,AgasIO)
        #fIOVecT[:,t+1] = ftotIO_tp1 #store the filtered parameters from previous iteration
        fIOVecT = hcat(fIOVecT, ftotIO_tp1 )

        ftotIO_t = ftotIO_tp1
        loglike += loglike_t
    end
    if dgp
        return Y_T , fIOVecT
    else
        return fIOVecT, loglike
    end
    end

score_driven_filter_or_dgp(Model::SdErgmDirBin1) = score_driven_filter_or_dgp(Model,array_2_vec_all_par(Model,Model.Par))
function gasScoreSeries( Model::SdErgmDirBin1,
                                dynPar_T::Array{<:Real,2};
                                obsT::Array{<:Real,2}=Model.obsT)
    N2,T = size(obsT);N = round(Int,N2/2)
    # Organize parameters of the GAS update equation
    sIO_T = ones(Real,2N,T)
    gradIO_T = ones(Real,2N,T)
    # identify(Model,UMallNodesIO)
    indTvNodesIO  = trues(2N)
    for t=1:T
        ftotIO_t = dynPar_T[:,t]
        degsIO_t = obsT[:,t] # vector of in and out degrees

             fIO_t = ftotIO_t[indTvNodesIO] #Time varying fitnesses
             thetas_mat_t_exp, exp_mat_t = StaticNets.expMatrix2(StaticNets.fooErgmDirBin1,ftotIO_t)
             exp_degIO_t = [sumSq(exp_mat_t,2);sumSq(exp_mat_t,1)]

             #The following part is a faster version of logLikelihood_t()
             loglike_t = sum(ftotIO_t.*degsIO_t)  .-   sum(log.(1 .+ thetas_mat_t_exp))
             # compute the score of observations at time t wrt parameters at tm1
             # THIS IS A NEW VERSION!!!!!!!!!!!!!!!!!!!!!!1------------------------------------------------------
             gradIO_t = sIO_t = identify(Model, (degsIO_t .- exp_degIO_t ) ) # only for TV paramaters
             # updating direction without rescaling (for the moment)
             I_t = scalingMatGas(Model,exp_mat_t)
             # display(diag(I_t))
             # display(size(gradIO_t))
             # display(size(exp_mat_t ))
              sIO_t = identify(Model, I_t\gradIO_t )
              sIO_t[gradIO_t.==0] .= 0
              sIO_T[:,t] =  sIO_t
              gradIO_T[:,t] =  gradIO_t
    end
    return sIO_T,gradIO_T
    end

sampl(Mod::SdErgmDirBin1,T::Int)=( N = length(Mod.groupsInds[1]) ;
                                        tmpTuple = score_driven_filter_or_dgp(Mod,[Mod.Par[1];Mod.Par[2];Mod.Par[3]];  dgpNT = (N,T) );
                                        degsIO_T = [sumSq(tmpTuple[1],3) sumSq(tmpTuple[1],2)];
                                        (SdErgmDirBin1(degsIO_T,Mod.Par,Mod.groupsInds,Mod.scoreScalingType),tmpTuple[1],tmpTuple[2]) )

# Estimation

function estSingSnap(Model::SdErgmDirBin1, degs_t::Array{<:Real,1}; groupsInds = Model.groupsInds, targetErr::Real=targetErrValDynNets)
    hatUnPar,~,~ = StaticNets.estimate(StaModType(Model); degIO = degs_t ,groupsInds = groupsInds[1], targetErr =  targetErr)
    return hatUnPar
 end
function estimateSnapSeq(Model::SdErgmDirBin1; degsIO_T::Array{<:Real,2}=Model.obsT,
                            targetErr::Real=1e-5,identPost=false,identIter=false)
    #this funciton does  not allow groups estimate a sequence of single node's
    # fitnesses
    return  StaticNets.estimate(StaticNets.SnapSeqNetDirBin1(degsIO_T),targetErr = targetErr ,identPost=identPost,identIter= identIter)
 end

function estimateTarg(Model::SdErgmDirBin1; SSest::Array{<:Real,2} =zeros(2,2),
             groupsInds::Array{<:Array{<:Real,1},1} = Model.groupsInds)
    "Estimate the GAS parameters of the GAS  model,
    using a targeting on the unconditional means. Then return the GASparameters "
    # Not ready for groups indices
    # only one B and one A
    N2,T = size(Model.obsT);N = round(Int,N2/2)
    println(N)
    NGW,GBA,ABgroupsIndNodesIO,indTvNodesIO = NumberOfGroupsAndABindNodes(Model,groupsInds)

    optims_opt, algo = setOptionsOptim(Model)

    #target the unconditional means from the sequence of single snapshots estimate
    if sum(SSest) == 0
        uncMeansIO =  meanSq(estimateSnapSeq(Model),2)
    else
        uncMeansIO =  meanSq(SSest,2)
    end

    uncMeansIO = identify(Model,uncMeansIO)
    ftotIO_0  = meanSq(SSest[:, 1:5],2)

    #functions needed if I want to enable grouping of nodes again
    # function WGroupsFromB(B::Array{<:Real,1})
    #      BGroupsW_IO = B[[groupsInds[2];groupsInds[2] ] ]
    #      WGroupsIO = uncMeansWGroupsIO.*(1-BGroupsW_IO)
    #      return WGroupsIO
    # end
    function restrAB(vecUnPar::Array{<:Real,1}) #restrict a vector of only A and B
         diag_B_Un = vecUnPar[1:GBA]
         diag_A_Un = vecUnPar[GBA+1:end]

         diag_B_Re = 1 ./ (1 .+ exp.( .- diag_B_Un))
         diag_A_Re = exp.(diag_A_Un)
         vecRePar =  [diag_B_Re; diag_A_Re]
         return vecRePar
         end
    function UnrestrAB(vecRePar::Array{<:Real,1}) #restrict a vector of only A and B
         diag_B_Re = vecRePar[1:GBA]
         diag_A_Re = vecRePar[GBA+1:end]
         diag_B_Un = log.(diag_B_Re ./ (1 .- diag_B_Re ))
         diag_A_Un = log.(diag_A_Re)
         vecUnPar =  [diag_B_Un; diag_A_Un]
         return vecUnPar
    end
    # #set the starting points for the optimizations
    B0_Re  = 0.99
    A0_Re_H  = 0.1
    A0_Re_L  = 0.000001
    if true && GBA>1
        #calcola autocorrelazioni dello score e usale per definire A_0
        meanDegsIO = meanSq(Model.obsT,2)
        constParsIO,~,~ = StaticNets.estimate(StaticNets.fooErgmDirBin1,degIO = meanDegsIO)
         fooPar = [zeros(2N), zeros(2N) ,zeros(2N)]
          indsGroups = [Int.(1:2N),Int.(1:2N)]
        modGasDirBin1 = DynNets.SdErgmDirBin1(Model.obsT, fooPar, indsGroups,"FISHER-DIAG")
        sIO_T,gIO_T =  gasScoreSeries(modGasDirBin1,repeat(constParsIO,1,T);obsT = Model.obsT)
        indsLikelyTV = dropdims(autocor(Float64.(gIO_T'),[1]).>0.2, dims = 1)
        Avec = ones(Real,GBA).*A0_Re_L
        println(size(indsLikelyTV))
        Avec[indsLikelyTV] .= A0_Re_H
    else
        Avec = ones(Real,GBA).*A0_Re_L
    end
    vGasParAll0_Un = UnrestrAB( [ones(Real,GBA).*B0_Re;Avec ]) # vGasParGroups0_Un#

    # objective function for the optimization
    function objfunGas(vecUnAB::Array{<:Real,1})# a function of the groups parameters
         ReB = 1 ./ (1 .+ exp.(.-vecUnAB[1:GBA]))
         ReA = exp.(vecUnAB[GBA+1:end])
         WNodes = uncMeansIO.*((maxLargeVal .+ ReB[1])/maxLargeVal)
        # println(typeof(WNodes))
         WNodes[indTvNodesIO] = uncMeansIO[indTvNodesIO].*(1 .- ReB[ABgroupsIndNodesIO[indTvNodesIO]])

         #StaticNets.expMatrix2(StaticNets.fooErgmDirBin1,WNodes)
         vecReGasPar = [WNodes;ReB;ReA ]
         #StaticNets.expMatrix2(StaticNets.fooErgmDirBin1,vecReGasPar[1:NGW])
         #NON POSSONO STARE NELLO STESSO VETTORE PARAMETRI DA OTTIMIZZARE DI TIPO FORWARDDIFF E TYPES NORMALI  ?      ?      ?      ?
         foo,loglikelValue = score_driven_filter_or_dgp( Model,vecReGasPar;
                                        groupsInds = groupsInds, ftotIO_0=ftotIO_0)
         return - loglikelValue
    end
    diff_mode = "forward"
    if diff_mode == "forward"
        if uppercase(Model.scoreScalingType) == "FISHER-EWMA"
            ADobjfunGas = objfunGas
        else
            ADobjfunGas = TwiceDifferentiable(objfunGas, vGasParAll0_Un; autodiff = :finite);
        end
    elseif diff_mode in ["reverse",  "backward"]

    end
    println(objfunGas(vGasParAll0_Un))
    optim_out2  = optimize(ADobjfunGas,vGasParAll0_Un,algo,optims_opt)
    outParAB = Optim.minimizer(optim_out2)
    outReAB = restrAB(outParAB)
    vGasParGroupsHat_Re = zeros(NGW + 2GBA)
    vGasParGroupsHat_Re[(1:NGW)[.!indTvNodesIO]] = uncMeansIO[.!indTvNodesIO]
    #println( uncMeansIO[.!indTvNodesIO])
    vGasParGroupsHat_Re[(1:NGW)[indTvNodesIO]] = uncMeansIO[indTvNodesIO].*(1 .- outReAB[1:GBA][ABgroupsIndNodesIO[indTvNodesIO]])
    vGasParGroupsHat_Re[1+NGW:end] = outReAB
    #println(vGasParGroupsHat_Re)
    arrayGasParHat_Re = vec_2_array_all_par(Model,vGasParGroupsHat_Re;groupsInds = groupsInds )
    conv_flag =  Optim.converged(optim_out2)
    # println(optim_out2)
    return  arrayGasParHat_Re, conv_flag
end

function estimateOld(Model::SdErgmDirBin1;start_values = [zeros(Real,10),zeros(Real,3),zeros(Real,3) ])
    T,N2 = size(Model.obsT);N = round(Int,N2/2)
    groupsInds = Model.groupsInds
    NGW,GBA,ABgroupsIndNodesIO,indTvNodesIO = NumberOfGroupsAndABindNodes(Model,groupsInds)

    optims_opt, algo = setOptionsOptim(Model)

    #set the starting points for the optimizations if not given
    if sum([sum(start_values[i]) for i =1:3 ])==0
         B0_ReGroups = 0.9 * ones(GBA)
         A0_ReGroups = 0.01 * ones(GBA)

         BgasNodes_0   = B0_ReGroups[ABgroupsIndNodesIO[indTvNodesIO]]
         AgasNodes_0   = A0_ReGroups[ABgroupsIndNodesIO[indTvNodesIO]]
         # choose the starting points for W in order to fix the Unconditional Means
         # equal to the static estimates for the average degrees(over time)
         MeanDegsT =  meanSq(Model.obsT,1)

         UMstaticGroups = estSingSnap(Model ,MeanDegsT; groupsInds = groupsInds )
         conv_flag1 = true

         GBA>0  ?    ( W0_Groups = UMstaticGroups.*(1 .- B0_ReGroups[[groupsInds[2];groupsInds[3]]])) : W0_Groups = UMstaticGroups
    else
         W0_Groups = start_values[1]
         B0_ReGroups = start_values[2]
         A0_ReGroups = start_values[3]
         conv_flag1 = true
    end
    if GBA == 0  return [UMstaticGroups, zeros(1), zeros(1)] ,  conv_flag1 end;
    # Estimate  jointly all  Gas parameters
    vGasParAll0_Un = unRestrictGasPar(Model,[W0_Groups ;B0_ReGroups;A0_ReGroups ] ) # vGasParGroups0_Un#
    # objective function for the optimization
    function objfunGas(vparUn::Array{<:Real,1})# a function of the groups parameters
              foo,loglikelValue = score_driven_filter_or_dgp( Model,restrictGasPar(Model,vparUn))
           return - loglikelValue
    end

    #Run the optimization without AD if a matrix sqrt is required in the filter
    if uppercase(Model.scoreScalingType) == "FISHER-EWMA"
        ADobjfunGas = objfunGas
    else
        ADobjfunGas = TwiceDifferentiable(objfunGas, vGasParAll0_Un; autodiff = :forward);
    end

    optim_out2  = optimize(ADobjfunGas,vGasParAll0_Un,algo,optims_opt)
    outParAll = Optim.minimizer(optim_out2)
    vGasParGroupsHat_Re = restrictGasPar(Model,outParAll)
    arrayGasParHat_Re = vec_2_array_all_par(Model,vGasParGroupsHat_Re )
    conv_flag = conv_flag1*Optim.converged(optim_out2)
    #println(optim_out2)
    return  arrayGasParHat_Re, conv_flag
    end

function estimateTargOld(Model::SdErgmDirBin1)
    "Estimate the GAS parameters of the GAS  model,
    using a targeting on the unconditional means. Then return the GASparameters "
    T,N2 = size(Model.obsT);N = round(Int,N2/2)
    groupsInds = Model.groupsInds
    NGW,GBA = NumberOfGroups(Model,groupsInds)
    optims_opt, algo = setOptionsOptim(Model)

    #target the unconditional mean
    ParTvSnap = estimateSnapSeq(Model) # sequence of snapshots estimate
    uncMeansWGroupsIO = dropdims(mean(ParTvSnap,1), dims =1)
    #functions needed for the specific targeting case
    function WGroupsFromB(B::Array{<:Real,1})
         BGroupsW_IO = B[[groupsInds[2];groupsInds[2] ] ]
         WGroupsIO = uncMeansWGroupsIO.*(1 .- BGroupsW_IO)
         return WGroupsIO
    end
    function restrAB(vecUnPar::Array{<:Real,1}) #restrict a vector of only A and B
         diag_B_Un = vecUnPar[1:GBA]
         diag_A_Un = vecUnPar[GBA+1:end]

         diag_B_Re = 1 ./ (1 .+ exp.(.-diag_B_Un))
         diag_A_Re = exp.(diag_A_Un)
         vecRePar =  [diag_B_Re; diag_A_Re]
         return vecRePar
         end
    function UnrestrAB(vecRePar::Array{<:Real,1}) #restrict a vector of only A and B
         diag_B_Re = vecRePar[1:GBA]
         diag_A_Re = vecRePar[GBA+1:end]
         diag_B_Un = log.(diag_B_Re ./ (1 .- diag_B_Re ))
         diag_A_Un = log.(diag_A_Re)
         vecUnPar =  [diag_B_Un; diag_A_Un]
         return vecUnPar
    end
    if GBA == 0  return [uncMeansWGroupsIO, zeros(1), zeros(1)] , Optim.converged(optim_out) end;
    # Estimate As and Bs

    # #set the starting points for the optimizations
    B0_ReGroups = 0.9 * ones(GBA)
    A0_ReGroups = 0.01 * ones(GBA)
    vGasParAll0_Un = UnrestrAB([B0_ReGroups;A0_ReGroups ] ) # vGasParGroups0_Un#
    # objective function for the optimization
    function objfunGas(vecUnAB::Array{<:Real,1})# a function of the groups parameters
         vecReAB = restrAB(vecUnAB)
         ReB = vecReAB[1:GBA]
         WGroups = WGroupsFromB(ReB)
         vecReGasPar = [WGroups;vecReAB]
         foo,loglikelValue = score_driven_filter_or_dgp( Model,vecReGasPar )
         return - loglikelValue
    end
    #Run the optimization
    if uppercase(Model.scoreScalingType) == "FISHER-EWMA"
        ADobjfunGas = objfunGas
    else
        ADobjfunGas = TwiceDifferentiable(objfunGas, vGasParAll0_Un; autodiff = :forward);
    end

    optim_out2  = optimize(ADobjfunGas,vGasParAll0_Un,algo,optims_opt)
    outParAB = Optim.minimizer(optim_out2)
    outReAB = restrAB(outParAB)
    vGasParGroupsHat_Re = [WGroupsFromB(outReAB[1:GBA]);outReAB]
    arrayGasParHat_Re = vec_2_array_all_par(Model,vGasParGroupsHat_Re )
    conv_flag =  Optim.converged(optim_out2)
    # println(optim_out2)
    return  arrayGasParHat_Re, conv_flag
    end

function forecastEvalGasNetDirBin1(obsNet_T::BitArray{3}, gasParEstOnTrain::Array{Array{Float64,1},1},Ttrain::Int;thVarConst = 0.0005,mat2rem=falses(obsNet_T[:,:,1]) )
    # this function evaluates the forecasting performances of the GAS model over the train sample.
    # to do so uses gas parameters previously estimated and observations at time t of the test sample
    # to forecast tv par at t+1, then using these estimates and t+1 observations goes to t+2 and so on

    @show N2 = length(gasParEstOnTrain[1]);N = round(Int,N2/2)
    T = length(obsNet_T[1,1,:])
    @show Ttest = T-Ttrain
    testObsNet_T = obsNet_T[:,:,Ttrain+1:end]
    degsIO_T = [sumSq(obsNet_T,2); sumSq(obsNet_T,1)]

    #disregard predictions for degrees that are constant in the training sample
    # gas shoud manage but AR fitness no and I dont want to advantage gas
    # to do so i need to count the number of links that are non constant (nnc)
    Nc,isConIn,isConOut = StaticNets.defineConstDegs(degsIO_T[:,1:Ttrain];thVarConst = thVarConst )
    noDiagIndnnc = putZeroDiag(((.!isConIn).*(.!isConOut')).* (.!mat2rem)    )

    @show Nlinksnnc =sum(noDiagIndnnc)
    # forecast fitnesses using Gas parameters and observations
    foreFit,~ = score_driven_filter_or_dgp( DynNets.SdErgmDirBin1(degsIO_T),[gasParEstOnTrain[1];gasParEstOnTrain[2];gasParEstOnTrain[3]])

    TRoc = Ttest-1
    #storage variables
    foreVals = 100 * ones(Nlinksnnc*(TRoc))
    realVals = trues( Nlinksnnc*(TRoc))

    lastInd=1

    for t=2:TRoc
        adjMat = testObsNet_T[:,:,t]   #[testObsNet_T[50:end,:,t];testObsNet_T[1:49,:,t]] #a shuffle to test that the code does not not work
        adjMat_tm1 = testObsNet_T[:,:,t-1]
        indZeroR = sum(adjMat_tm1,dims = 2).==0
        indZeroC = sum(adjMat_tm1,dims = 1).==0
        indZeroMat = indZeroR.*indZeroC
    #    println(sum(indZeroMat)/(length(indZeroC)^2))
        expMat = StaticNets.expMatrix(StaticNets.fooErgmDirBin1,foreFit[:,Ttrain+1:end][:,t])

        #expMat[indZeroMat] = 0
        foreVals[lastInd:lastInd+Nlinksnnc-1] = expMat[noDiagIndnnc]
        realVals[lastInd:lastInd+Nlinksnnc-1] = adjMat[noDiagIndnnc]
        lastInd += Nlinksnnc
    end
    return realVals,foreVals,foreFit
 end

function forecastEvalGasNetDirBin1(obsMat_allT::BitArray{3}, Ttrain::Int;thVarConst = 0.0005,mat2rem=falses(obsNet_T[:,:,1]) )
    # if all observations and the length of training sample are given as input, then
    # estimate the gas par and run the one step forecast
    degsIO_T = [sumSq(obsMat_allT,2);sumSq(obsMat_allT,1)]
    N2,T = size(degsIO_T);N = round(Int,N2/2)
    Ttest = T - Ttrain
    #Define the train model and estimate
    modGasDirBin1_eMidTrain = DynNets.SdErgmDirBin1(degsIO_T[:,1:Ttrain])
    estTargDirBin1_eMidTrain,~ = DynNets.estimateTarg(modGasDirBin1_eMidTrain)
    gasParEstOnTrain = estTargDirBin1_eMidTrain
    #produce a vector of real values and probabilities from 1 step ahead forecasting
    foreEval = forecastEvalGasNetDirBin1(obsMat_allT,gasParEstOnTrain,Ttrain,thVarConst = thVarConst,mat2rem =mat2rem)
    tp,fp,AUC = Utilities.rocCurve(foreEval[1],foreEval[2])
    return tp,fp,AUC,foreEval
 end

function multiSteps_predict_score_driven_par( Model::SdErgmDirBin1,N::Int,degsIO_t_in::Array{<:Real,1},
                         ftotIO_t0_input::Array{<:Real,1},
                         WgasIO_in::Array{<:Real,1},BgasIO_in::Array{<:Real,1},AgasIO_in::Array{<:Real,1},Nsample::Int,Nsteps::Int )
     #Simula Nsteps nel futuro per Nsample volte
     degsIO_t0 = copy(degsIO_t_in)
     ftotIO_t0 = copy(ftotIO_t0_input)
     WgasIO = copy(WgasIO_in)
     BgasIO = copy(BgasIO_in)
     AgasIO = copy(AgasIO_in)
     ftotIO_T = zeros(length(ftotIO_t0),Nsteps)

     indTvNodesIO = Model.groupsInds[2].!=0
     expMat_T_means = zeros(N,N,Nsteps)
     expMat_T_Y = zeros(N,N,Nsteps)
     ftotIO_t1  ,~ = predict_score_driven_par(Model,N,degsIO_t0,ftotIO_t0,eye(N,N),
            indTvNodesIO,WgasIO,BgasIO,AgasIO)
     ftotIO_T[:,1] = copy(ftotIO_t1)
     for i=1:Nsample
         ftotIO_t1_tmp = copy(ftotIO_t1)
         expMat_t = StaticNets.expMatrix(StaModType(Model),ftotIO_t1_tmp )
         expMat_T_means[:,:,1] = expMat_t
         Y_t = StaticNets.samplSingMatCan(StaModType(Model),expMat_t)
         degsIO_t = [dropdims(sum(Y_t,dims = 2),dims =2); dropdims(sum(Y_t,dims = 1), dims =1)]
         #matrice media dell matrici campionate
         expMat_T_Y[:,:,1] = expMat_T_Y[:,:,1]  .+  Y_t ./ Nsample
         #il primo step viene fatto usando l'osservzione reale
         ftotIO_tp1  ,~ = predict_score_driven_par(Model,N,degsIO_t,ftotIO_t1_tmp,eye(N,N),
                            indTvNodesIO,WgasIO,BgasIO,AgasIO)
         ftotIO_t = copy(ftotIO_tp1)
         ftotIO_T[:,1] = copy(ftotIO_t)
         expMat_t = StaticNets.expMatrix(StaModType(Model),ftotIO_t )
        # println(sum(ftotIO_t0))

         for t=2:Nsteps
            # fai update dei parametri
             ftotIO_tp1  ,~ = predict_score_driven_par(Model,N,degsIO_t,ftotIO_t,eye(expMat_t),
                                                 indTvNodesIO,WgasIO,BgasIO,AgasIO)
             ftotIO_t = copy(ftotIO_tp1)
             # da ora stai usando i parametri  t
             ftotIO_T[:,t] = ftotIO_T[:,t] .+ ftotIO_t ./ Nsample

             expMat_t = StaticNets.expMatrix(StaModType(Model),ftotIO_t )
            # t==2  ?     println(mean(expMat_t)) : ()
             #matrice media delle matrici medie
             expMat_T_means[:,:,t] = expMat_T_means[:,:,t]  .+  expMat_t ./ Nsample

             Y_t = StaticNets.samplSingMatCan(StaModType(Model),expMat_t)
             #matrice media dell matrici campionate
             expMat_T_Y[:,:,t] = expMat_T_Y[:,:,t]  .+  Y_t ./ Nsample
             degsIO_t = [dropdims(sum(Y_t,dims = 2), dims =2); dropdims(sum(Y_t,dims = 1), dims =1)]

         end

     end
     return expMat_T_means, ftotIO_T,expMat_T_Y
  end

function multiStepsForecastExpMat( Model::SdErgmDirBin1,obsNet_T::BitArray{3},
      gasParEstOnTrain::Array{Array{Float64,1},1},Nsample::Int,Nsteps::Int,Ttrain::Int;meanOfMeans=false)

      N2 = length(gasParEstOnTrain[1]);N = round(Int,N2/2)
      T = length(obsNet_T[1,1,:])
      degsIO_T = [sumSq(obsNet_T,2);sumSq(obsNet_T,1)]
      foreFitGas1,~ = DynNets.score_driven_filter_or_dgp( Model,
      [gasParEstOnTrain[1];gasParEstOnTrain[2];gasParEstOnTrain[3]])
      Ttest = T-Ttrain

      expMat_T_means = zeros(N,N,T)
      expMat_T_Y = zeros(N,N,T)
      foreFit_T = zeros(2N,T)
      expMat_T_meanPar = zeros(N,N,T)
      for t=Ttrain +1:T
          tFit = t-Nsteps
          tmpAllMat_means,tmpFit,tmpAllMat_Y = multiSteps_predict_score_driven_par(Model,N,degsIO_T[:,tFit],foreFitGas1[:,tFit],
                                gasParEstOnTrain[1],gasParEstOnTrain[2],gasParEstOnTrain[3],Nsample,Nsteps)
          expMat_T_means[:,:,tFit+Nsteps] = tmpAllMat_means[:,:,Nsteps]
          expMat_T_Y[:,:,tFit+Nsteps] = tmpAllMat_Y[:,:,Nsteps]
          foreFit_T[:,tFit+Nsteps] = tmpFit[:,end]
           expMat_T_meanPar[:,:,tFit+Nsteps] =   StaticNets.expMatrix(StaModType(Model),tmpFit[:,end])
      end
      return expMat_T_means,foreFit_T,expMat_T_Y,expMat_T_meanPar
  end

function multiStepsForecastExpMat_roll( Model::SdErgmDirBin1,obsNet_T::BitArray{3},
      gasEst_rolling_input::Array{Array{Array{Float64,1},1},1},gasFiltAndForeFitFromRollEst_input::Array{Float64,2},Nsample::Int,Nsteps::Int,Ttrain::Int;meanOfMeans=false)

      gasEst_rolling = copy(gasEst_rolling_input)
      gasFiltAndForeFitFromRollEst = copy(gasFiltAndForeFitFromRollEst_input)
      N2 = length(gasEst_rolling[1][1]);N = round(Int,N2/2)
      T = length(obsNet_T[1,1,:])
      degsIO_T = [sumSq(obsNet_T,2);sumSq(obsNet_T,1)]
      Ttest = T-Ttrain

      expMat_T_means = zeros(N,N,T)
      expMat_T_Y = zeros(N,N,T)
      foreFit_T = zeros(2N,T)
      expMat_T_meanPar = zeros(N,N,T)
      # forecast the parameters at each t between Ttrain+1 and T using the observation at t-N_Nsteps
      # and the static gas estimates obtained from  t-Nsteps-Ttrain:T-Nsteps
      for t=Ttrain +1:T
          tEst = max(t-Nsteps - Ttrain,1)
          tObs = t-Nsteps-1

          tmpAllMat_means,tmpFit,tmpAllMat_Y = multiSteps_predict_score_driven_par(Model,N,degsIO_T[:,tObs],gasFiltAndForeFitFromRollEst[:,tObs+1],
                                gasEst_rolling[tEst][1],gasEst_rolling[tEst][2],gasEst_rolling[tEst][3],Nsample,Nsteps)
          expMat_T_means[:,:,tObs+Nsteps] = tmpAllMat_means[:,:,Nsteps]
          expMat_T_Y[:,:,tObs+Nsteps] = tmpAllMat_Y[:,:,Nsteps]
          foreFit_T[:,tObs+Nsteps] = tmpFit[:,end]
           expMat_T_meanPar[:,:,tObs+Nsteps] =   StaticNets.expMatrix(StaModType(Model),tmpFit[:,end])
      end
      return expMat_T_means,foreFit_T,expMat_T_Y,expMat_T_meanPar
  end

function logLike_t(Model::SdErgmDirBin1, obsT, vReGasPar)
 groupsInds = Model.groupsInds
      N2,T = size(obsT);N = round(Int,N2/2)
      NGW,GBA,ABgroupsIndNodesIO,indTvNodesIO = NumberOfGroupsAndABindNodes(Model, groupsInds)
      # Organize parameters of the GAS update equation
      WGroupsIO = vReGasPar[1:NGW]
      #    StaticNets.expMatrix2(StaticNets.fooErgmDirBin1,WGroupsIO )
      W_allIO = WGroupsIO[groupsInds[1]]
      BgasGroups  = vReGasPar[NGW+1:NGW+GBA]
      AgasGroups  = vReGasPar[NGW+GBA+1:NGW+2GBA]
      #distribute nodes among groups with equal parameters
      AgasIO = AgasGroups[ABgroupsIndNodesIO[indTvNodesIO]]
      BgasIO = BgasGroups[ABgroupsIndNodesIO[indTvNodesIO]]
      WgasIO   = W_allIO[indTvNodesIO]
      UMallNodesIO = W_allIO
      UMallNodesIO  =  WgasIO ./ (1 .- BgasIO)
      ftotIO_t =  identify(Model,UMallNodesIO)
      I_tm1 = UniformScaling(N)
      loglike_t = zero(Real)
      for t=1:T-1
            degsIO_t = obsT[:,t] # vector of in and out degrees
            ftotIO_tp1,loglike_t = predict_score_driven_par(Model,N,degsIO_t,ftotIO_t,I_tm1,
                                                  indTvNodesIO,WgasIO,BgasIO,AgasIO)
            ftotIO_t = ftotIO_tp1
      end
      return  loglike_t::T where T <:Real
end

function logLike_T(Model::SdErgmDirBin1, obsT, vReGasPar)
      groupsInds = Model.groupsInds
      N2,T = size(obsT);N = round(Int,N2/2)
      NGW,GBA,ABgroupsIndNodesIO,indTvNodesIO = NumberOfGroupsAndABindNodes(Model, groupsInds)
      # Organize parameters of the GAS update equation
      WGroupsIO = vReGasPar[1:NGW]
      #    StaticNets.expMatrix2(StaticNets.fooErgmDirBin1,WGroupsIO )
      W_allIO = WGroupsIO[groupsInds[1]]
      BgasGroups  = vReGasPar[NGW+1:NGW+GBA]
      AgasGroups  = vReGasPar[NGW+GBA+1:NGW+2GBA]
      #distribute nodes among groups with equal parameters
      AgasIO = AgasGroups[ABgroupsIndNodesIO[indTvNodesIO]]
      BgasIO = BgasGroups[ABgroupsIndNodesIO[indTvNodesIO]]
      WgasIO   = W_allIO[indTvNodesIO]
      # start values equal the unconditional mean, but  constant ones remain equal to the unconditional mean, hence initialize as:
      UMallNodesIO = W_allIO
      UMallNodesIO  =  WgasIO ./ (1 .- BgasIO)
      ftotIO_t =  identify(Model,UMallNodesIO)
      I_tm1 = UniformScaling(N)
      loglike_T = zero(Real)
      for t=1:T-1
            degsIO_t = obsT[:,t] # vector of in and out degrees
            ftotIO_tp1,loglike_t = predict_score_driven_par(Model,N,degsIO_t,ftotIO_t,I_tm1,
                                                  indTvNodesIO,WgasIO,BgasIO,AgasIO)
            ftotIO_t = ftotIO_tp1
            loglike_T += loglike_t
      end
      return  loglike_T::T where T <:Real
end


text:
# Illustrations for SolarIrradiance documentation
# Copyright (c) 2021 Pierre Haessig

using SolarIrradiance
using PyPlot
#pygui(true) # for interactive plot windows

## Load GHI data ##

# 2012, day 136: sun & clouds
year = 2012
n = 136
GHI_day = [0.0, 0.0, 0.0, 0.0, 0.0, 39.18, 209.0, 374.0, 240.0, 476.0, 379.0, 619.0, 872.01, 664.0, 680.0, 434.0, 311.0, 313.0, 169.0, 20.88, 0.0, 0.0, 0.0, 0.0]

# 2012, day 137: sun
n = 137
GHI_day = [0.0, 0.0, 0.0, 0.0, 0.0, 64.0, 218.0, 396.0, 566.0, 709.0, 827.01, 892.01, 851.01, 834.01, 756.01, 666.0, 519.0, 346.0, 168.0, 21.88, 0.0, 0.0, 0.0, 0.0]

## Custom colormaps for sets of license

# https://gka.github.io/palettes/#/9|s|0030dd,ef272d|ffffe0,ff005e,93003a|1|1
colorsBR = ["#0030dd", "#622cc7", "#8628b1", "#9f259b", "#b32486", "#c52370", "#d4235b", "#e22545", "#ef272d"]
cmapBR = matplotlib.colors.LinearSegmentedColormap.from_list("BlueRed", colorsBR)

# https://gka.github.io/palettes/#/9|s|3b2fa7,007300,887600,eb212a|ffffe0,ff005e,93003a|0|0
colorsBGR = ["#3b2fa7", "#254968", "#0f622a", "#117300", "#447400", "#777600", "#a1610b", "#c6411a", "#eb212a"]
cmapBGR = matplotlib.colors.LinearSegmentedColormap.from_list("BlueGreenRed", colorsBGR)

# https://gka.github.io/palettes/#/9|d|0030dd,4d4d4d|4d4d4d,df0822|0|0
colorsBKR = ["#0030dd", "#1337b9", "#263e95", "#3a4671", "#4d4d4d", "#713c42", "#962a37", "#ba192d", "#df0822"]
# https://gka.github.io/palettes/#/9|d|0024cf,505050|505050,df0822|0|1
colorsBKR = ["#0024cf", "#3d31ae", "#4d3d8e", "#52476f", "#505050", "#794c45", "#9c443a", "#be342e", "#df0822"]
cmapBKR = matplotlib.colors.LinearSegmentedColormap.from_list("BlueGrayRed", colorsBKR)

# https://gka.github.io/palettes/#/9|d|00aa00,b60053|b60053,0000aa|0|1
colorsGRB = ["#00aa00", "#659129", "#89743b", "#a25148", "#b60053", "#9e0069", "#82007e", "#5e0094", "#0000aa"]
cmapGRB = matplotlib.colors.LinearSegmentedColormap.from_list("GreenRedBlue", colorsGRB)

## Transpose GHI for a given slope & azimuth ##

# Time of the measurement
tc_range = 0:23
dt = 1/60 # 1 minute ~ instantaneous (0.0 is orbidden)

# Location (Rennes, France):
lat = 48.117 #°
lon = -1.678 #°

# Panel orientation:
slope = 40 #°
azimuth = 90 #° → west-facing
albedo = 0.5; # in [0,1]

Gplane = zeros(length(tc_range))

for (k, tck) in enumerate(tc_range)
    GHIk = GHI_day[k]
    Gplane[k] = global_radiation_tilt(GHIk, n, tck, dt, lat, lon, slope, azimuth, albedo)
end

## Plot the single transpose output ##

fig, ax = subplots(figsize=(5,3.5))
ax.plot(tc_range, GHI_day, "o-", label="horizontal (data)")
ax.plot(tc_range, Gplane, "D-", label="tilted (estimate)")
ax.grid(true)
ax.set(
    title="Irradiance transposition\non a west-facing plane tilted by $(slope)°",
    xlabel="time (UTC hours)",
    ylabel="Irradiance (W/m²)",
)
ax.legend()
fig.tight_layout()
fname = "GHI_$(year)-$(n)_trans_s$(slope)_a$(azimuth).png"
fig.savefig(fname, dpi=300)

## Plot a range of azimuths ##

azi_range = -90:15:90

Gplane_azi = zeros(length(tc_range), length(azi_range))

for (i, azimuth) in enumerate(azi_range)
    for (k, tck) in enumerate(tc_range)
        GHIk = GHI_day[k]
        Gplane_azi[k,i] = global_radiation_tilt(GHIk, n, tck, dt, lat, lon, slope, azimuth, albedo)
    end
end

# Plot azimuth range
fig, ax = subplots(figsize=(5,3.5))

ax.plot(tc_range, GHI_day, "-", label="GHI data",
        color=(0,0.7,0, 0.5), ms=3, lw=3, alpha=0.5)

for (i, azimuth) in enumerate(azi_range)
    lw=1
    ls="-"
    if azimuth == -90
        label = "−90° (E)"
    elseif azimuth == 0
        label = "0° (S)"
    elseif azimuth == +90
        label = "+90° (W)"
    else
        label = ""
        lw=0.5
        ls="-"
    end
    color = cmapBKR((i-1)/(length(azi_range)-1))
    ax.plot(tc_range, Gplane_azi[:,i], ls, ms=3,
            color = color, label=label, lw=lw)
end

ax.grid(true)
ax.set(
    title="Irradiance transposition\non a plane tilted by $(slope)° with varying azimuths",
    xlabel="time (UTC hours)",
    ylabel="Irradiance (W/m²)",
)
ax.set_ylim(0, 1400)
ax.legend(ncol=2)
fig.tight_layout()
fname = "GHI_$(year)-$(n)_trans_s$(slope)_azim_range.png"
fig.savefig(fname, dpi=300)

## Plot a range of slopes ##

azimuth = 0 #°
slo_range = 0:15:90

Gplane_slo = zeros(length(tc_range), length(slo_range))

for (i, slope) in enumerate(slo_range)
    for (k, tck) in enumerate(tc_range)
        GHIk = GHI_day[k]
        Gplane_slo[k,i] = global_radiation_tilt(GHIk, n, tck, dt, lat, lon, slope, azimuth, albedo)
    end
end

# Plot slope range
fig, ax = subplots(figsize=(5,3.5))

ax.plot(tc_range, GHI_day, "-", label="GHI data",
        color=(0,0.7,0, 0.5), ms=3, lw=3, alpha=0.5)

for (i, slope) in enumerate(slo_range)
    lw=1
    ls="-"
    if slope == 0
        label = "0°"
    elseif slope == +45
        label = "45°"
    elseif slope == +90
        label = "90°"
    else
        label = ""
        lw=0.5
        ls="-"
    end
    color = cmapGRB((i-1)/(length(slo_range)-1))
    ax.plot(tc_range, Gplane_slo[:,i], ls, ms=3,
            color = color, label=label, lw=lw)
end

ax.grid(true)
ax.set(
    title="Irradiance transposition\non a south-facing plane with varying slopes",
    xlabel="time (UTC hours)",
    ylabel="Irradiance (W/m²)",
)
ax.set_ylim(0, 1400)
ax.legend(ncol=2, loc="upper center")
fig.tight_layout()
fname = "GHI_$(year)-$(n)_trans_a$(azimuth)_slope_range.png"
fig.savefig(fname, dpi=300)

## Animation of a range of azimuths ##

slope = 40 #°
azi_range = -90:3:90 # → 2s @ 30 fps

lineG_data = zeros(2, length(tc_range))
lineG_data[1,:] = tc_range

function update_Gplane_azi(i, lineG)
    azimuth = azi_range[i]
    # Compute Gplane
    for (k, tck) in enumerate(tc_range)
        GHIk = GHI_day[k]
        lineG_data[2,k] = global_radiation_tilt(GHIk, n, tck, dt, lat, lon, slope, azimuth, albedo)
    end
    ax.set_title("Irradiance transposition on a plane\n tilted by $(slope)°, azimuth $(azimuth)°")
    # Update plot
    lineG.set_data(lineG_data)
    return (lineG, )
end

# Figure creation with null Gplane
fig, ax = subplots(figsize=(5,3.5))
ax.plot(tc_range, GHI_day, "o-", label="horizontal (data)")
lineG = ax.plot(tc_range, lineG_data[2,:], "D-", label="tilted (estimate)")[1]
ax.grid(true)
ax.set(
    title="Irradiance transposition on a plane\n tilted by $(slope)°, azimuth XX°",
    xlabel="time (UTC hours)",
    ylabel="Irradiance (W/m²)",
    ylim=(0,1300)
)
ax.legend(ncol=2, loc="upper center")
fig.tight_layout()

# Creating the Animation object
azim_ani = matplotlib.animation.FuncAnimation(
    fig, update_Gplane_azi, 1:length(azi_range), fargs=(lineG,), interval=33.33)

fname = "GHI_$(year)-$(n)_trans_s$(slope)_azim_range.gif"
writer = matplotlib.animation.PillowWriter(fps=30)
fname = "GHI_$(year)-$(n)_trans_s$(slope)_azim_range.mp4"
writer = matplotlib.animation.FFMpegWriter(fps=30, metadata=Dict(
        "title" => "GHI transposition on a plane tilted by $(slope)°, for a range of azimuths",
        "artist" => "Pierre Haessig",
        "date" => "2021"))
azim_ani.save(fname, dpi=216, writer=writer) # 5"*216 = 1080 px

## Show azimuth of the plane

fig = plt.figure()
ax = fig.add_subplot(1,1,1)

function aziplot_setup(ax, plane_color = "tab:red", ra = 1.3)
    # N,E,S,W arrow
    NESW_color = "gray"
    ax.arrow(0, -ra , 0, 2ra, head_width=0.1, length_includes_head=true, color=NESW_color)
    ax.arrow(-ra , 0, 2ra, 0, head_width=0, color=NESW_color)
    ax.text(0, 1.1ra, "S\n0°", ha="center", va="bottom")
    ax.text(1.1ra, 0, "90° W", ha="left", va="center")
    ax.text(0, -1.1ra, "180°\nN", ha="center", va="top")
    ax.text(-1.1ra, 0, "E −90°", ha="right", va="center")

    # Plane orientation plot
    azimuth = 0 #°
    cosγ = cos(deg2rad(azimuth))
    sinγ = sin(deg2rad(azimuth))
    # plane
    plane_line = ax.plot([-cosγ, cosγ], [sinγ, -sinγ], color=plane_color, lw=3)[1]
    # normal arrow
    plane_arrow = ax.arrow(0, 0, sinγ, cosγ, head_width=0.1, length_includes_head=true, color=plane_color)
    plane_arrow_xy0 = plane_arrow.get_xy()

    # Axes setup:
    lim = 1.5ra
    ax.set(
        xlim=(-lim, lim),
        ylim=(-lim, lim),
        aspect = "equal"
    )
    ax.axis("off")

    return plane_arrow, plane_line, plane_arrow_xy0
end

"move plane to azimuth[i]"
function aziplot_update(i, plane_arrow, plane_line, plane_arrow_xy0)
    azimuth = azi_range[i]
    cosγ = cos(deg2rad(azimuth))
    sinγ = sin(deg2rad(azimuth))
    # plane
    plane_line.set_data([-cosγ, cosγ], [sinγ, -sinγ])
    # normal arrow
    rot = matplotlib.transforms.Affine2D().rotate_deg(-azimuth)
    xyr = rot.transform(plane_arrow_xy0)
    plane_arrow.set_xy(xyr)
    return plane_arrow, plane_line
end

plane_arrow, plane_line, plane_arrow_xy0 = aziplot_setup(ax)
aziplot_update(2, plane_arrow, plane_line, plane_arrow_xy0)

## Animation of a range of azimuth + azi_plot ##

slope = 40 #°
fps = 30
azi_range = -90:3:90 # → 2s @ 30 fps
print("azimuth anim $(azi_range[1])° to $(azi_range[end])°: $(length(azi_range)/fps)s @ $(fps) fps\n")

lineG_data = zeros(2, length(tc_range))
lineG_data[1,:] = tc_range

function update_Gplane_azi(i, lineG, ax)
    azimuth = azi_range[i]
    # Compute Gplane
    for (k, tck) in enumerate(tc_range)
        GHIk = GHI_day[k]
        lineG_data[2,k] = global_radiation_tilt(GHIk, n, tck, dt, lat, lon, slope, azimuth, albedo)
    end
    ax.set_title("Irradiance transposition on a plane\n tilted by $(slope)°, azimuth $(azimuth)°")
    # Update plot
    lineG.set_data(lineG_data)
    return (lineG, )
end

function update_azianim(i, lineG, axG, plane_arrow, plane_line, plane_arrow_xy0)
    lineG = update_Gplane_azi(i, lineG, axG)[1]
    plane_arrow, plane_line = aziplot_update(i, plane_arrow, plane_line, plane_arrow_xy0)
    return (lineG, plane_arrow, plane_line)
end

# Figure creation with null Gplane
fig = figure(figsize=(9, 4.5))

gs = fig.add_gridspec(1, 2, width_ratios=[3, 1])
ax1 = fig.add_subplot(gs[1])
ax2 = fig.add_subplot(gs[2])

plane_arrow, plane_line = aziplot_setup(ax2)

# GHI plot setup
ax1.plot(tc_range, GHI_day, "o-", label="horizontal (data)")
lineG = ax1.plot(tc_range, lineG_data[2,:], "D-", label="tilted (estimate)")[1]
ax1.grid(true)
ax1.set(
    title="Irradiance transposition on a plane\n tilted by $(slope)°, azimuth XX°",
    xlabel="time (UTC hours)",
    ylabel="Irradiance (W/m²)",
    ylim=(0,1300)
)
ax1.legend(ncol=2, loc="upper center")
fig.tight_layout()

# Creating the Animation object
azim_ani = matplotlib.animation.FuncAnimation(
    fig, update_azianim, 1:length(azi_range),
    fargs=(lineG, ax1, plane_arrow, plane_line, plane_arrow_xy0),
    interval=33.33)

# fname = "GHI_$(year)-$(n)_trans_s$(slope)_azim_range2.gif"
# writer = matplotlib.animation.PillowWriter(fps=30)
fname = "GHI_$(year)-$(n)_trans_s$(slope)_azim_range2.mp4"
writer = matplotlib.animation.FFMpegWriter(fps=30, metadata=Dict(
        "title" => "GHI transposition on a plane tilted by $(slope)°, for a range of azimuths",
        "artist" => "Pierre Haessig",
        "date" => "2021"))
azim_ani.save(fname, dpi=240, writer=writer) # 4.5"*240 = 1080 px

## Animation of a range of slopes ##

azimuth = 0 #°
slo_range = 0:1.5:90 # → 2s @ 30 fps

lineG_data = zeros(2, length(tc_range))
lineG_data[1,:] = tc_range

function update_Gplane_slo(i, lineG)
    slope = slo_range[i]
    # Compute Gplane
    for (k, tck) in enumerate(tc_range)
        GHIk = GHI_day[k]
        lineG_data[2,k] = global_radiation_tilt(GHIk, n, tck, dt, lat, lon, slope, azimuth, albedo)
    end
    ax.set_title("Irradiance transposition on a plane with azimuth $(azimuth)°,\n tilted by $(slope)°")
    # Update plot
    lineG.set_data(lineG_data)
    return (lineG, )
end

# Figure creation with null Gplane
fig, ax = subplots(figsize=(5,3.5))
ax.plot(tc_range, GHI_day, "o-", label="horizontal (data)")
lineG = ax.plot(tc_range, lineG_data[2,:], "D-", label="tilted (estimate)")[1]
ax.grid(true)
ax.set(
    title="Irradiance transposition on a plane with azimuth $(azimuth)°,\n tilted by XX°",
    xlabel="time (UTC hours)",
    ylabel="Irradiance (W/m²)",
    ylim=(0,1300)
)
ax.legend(ncol=2, loc="upper center")
fig.tight_layout()

# Creating the Animation object
azim_slo = matplotlib.animation.FuncAnimation(
    fig, update_Gplane_slo, 1:length(slo_range), fargs=(lineG,), interval=33.33)

fname = "GHI_$(year)-$(n)_trans_a$(azimuth)_slope_range.mp4"
azim_slo.save(fname, dpi=216, # 5"*216 = 1080 px
    metadata=Dict(
        "title" => "GHI transposition on a plane with azimuth $(azimuth)°, for a range of slopes",
        "artist" => "Pierre Haessig",
        "date" => "2021")
)

ext: jl


text:
using InterpretableValidation
using Distributions

x = collect(1:10.)
mvts = MvTimeseriesDistribution(:x => IID(x, Uniform()), :y => IID(x, Normal()), :z => IID(x, Categorical(5)))

function target_expr(data::Dict{Symbol, Array{Float64}})
    truth_vals = ((data[:x] .> 0.75) .& (data[:y] .< -0.75)) .| ((data[:x] .< 0.25) .& (data[:y] .> -0.25))
    truth_vals[1] = truth_vals[1] &&  (data[:z][1] == 5)
    return (length(truth_vals) - sum(truth_vals)) / length(truth_vals)
end

results = optimize(target_expr, mvts)
println("loss: ", results.loss, " expressions: ", results.expr)


text:
module ActivateDreamBlocks

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Activate Space Jams (Everest)" => Ahorn.EntityPlacement(
        Maple.ActivateDreamBlocksTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
# Regexes prefixed with 'r'

#=
+-------------------------+---------+----------------+
| MATCH TERM              | REGEX   | EXAMPLE        |
+-------------------------+---------+----------------+
| Any Char except \n      | .       | r"."           |
+-------------------------+---------+----------------+
| Character 'x'           | x       | r"x"           |
+-------------------------+---------+----------------+
| String "xy"             | xy      | r"xy"          |
+-------------------------+---------+----------------+
| x or y                  | x|y     | r"x|y"         |
+-------------------------+---------+----------------+
| Escape                  | \       | r"\n"          |
+-------------------------+---------+----------------+
| One char of a,b,c,d     | [ab-d]  | r"[ab-d]"      |
+-------------------------+---------+----------------+
| One char except a,b,c,d | [^ab-d] | r"[^ab-d]"     |
+-------------------------+---------+----------------+
| One digit               | \d      | r"\d"          |
+-------------------------+---------+----------------+
| One non digit           | \D      | r"\D"          |
+-------------------------+---------+----------------+
| one whitespace          | \s      | r"\s"          |
+-------------------------+---------+----------------+
| one non whitespace      | \S      | r"\S"          |
+-------------------------+---------+----------------+
| Newline                 | \n      | r"\n"          |
+-------------------------+---------+----------------+
| Tab                     | \t      | r"\t"          |
+-------------------------+---------+----------------+
| Null character          | \0      | r"\0"          |
+-------------------------+---------+----------------+
| 0+ occurrences          | *       | r"x*"          |
+-------------------------+---------+----------------+
| 1+ occurrences          | +       | r"x+"          |
+-------------------------+---------+----------------+
| 0-1 occurrences         | ?       | r"x?"          |
+-------------------------+---------+----------------+
| exactly 2 times         | {2}     | r"\d{2}"       |
+-------------------------+---------+----------------+
| between 2 and 4 times   | {2,4}   | r"x{2,4}"      |
+-------------------------+---------+----------------+
| 2+ times                | {2,}    | r"x{2,}"       |
+-------------------------+---------+----------------+
| Start of string         | ^       | r"^abc$"       |
+-------------------------+---------+----------------+
| End of string           | $       | r"abc$"        |
+-------------------------+---------+----------------+
| Word boundary           | \b      | r"abc\bdef\b." |
+-------------------------+---------+----------------+
| Not a word boundary     | \B      | r"\Babc\B"     |
+-------------------------+---------+----------------+
=#

regex = r"^\s*(?:#|$)";
println(typeof(regex));

# Regex check
println(occursin(regex, "#abcdef"));

# MATCH
println(match(regex, "#abcdef"));

# RegexMatch object info
m = match(regex, "#abcdef");

# substring matched
println(m.match);
# captured substring as an array
println(m.captures);
# offset where the match begins
println(m.offset);
# all offsets if multiple matches occur
println(m.offsets);

# Example
m = match(r"(?<hour>\d+):(?<minute>\d+)", "12:45");

println(m[1]);
println(m[:minute]);

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
# Misfit functions

## Predicates ##

@doc """`isstable(sys)`

Returns `true` if `sys` is stable, else returns `false`.""" ->
function isstable(sys::LTISystem)
    if sys.Ts == 0
        if any(real(pole(sys)).>=0)
            return false
        end
    else
        if any(abs(pole(sys)).>=1)
            return false
        end
    end
    return true
end

@doc """`iscontinuous(sys)`

Returns `true` if `sys` is continuous, else returns `false`.""" ->
function iscontinuous(sys::LTISystem)
    return sys.Ts == 0
end

@doc """`issiso(sys)`

Returns `true` if `sys` is SISO, else returns `false`.""" ->
function issiso(sys::LTISystem)
    return sys.ny == 1 && sys.nu == 1
end

@doc """`isproper(tf)`

Returns `true` if the `TransferFunction` is proper. This means that order(den)
\>= order(num))""" ->
function isproper(t::TransferFunction)
    for s in t.matrix
        if length(num(s)) > length(den(s))
            return false
        end
    end
    return true
end

## Helpers (not public) ##

# Convert the argument to a Matrix{Float64}
function float64mat(A::Vector)
    A = reshape(A, size(A, 1), 1)
    return float64mat(A)
end
float64mat(A::Matrix) = map(Float64,A)
float64mat(A::Matrix{Float64}) = A
float64mat(A::Real) = float64mat([A])

# Ensures the metadata for an LTISystem is valid
function validate_names(kwargs, key, n)
    names = get(kwargs, key, "")
    if names == ""
        return fill(String(""), n)
    elseif isa(names, Vector) && eltype(names) <: AbstractString
        return String[names[i] for i = 1:n]
    elseif isa(names, AbstractString)
        return String[names * "$i" for i = 1:n]
    else
        error("$key must be of type `AbstractString` or Vector{AbstractString}")
    end
end

# Format the metadata for printing
function format_names(names::Vector{String}, default::AbstractString, unknown::AbstractString)
    n = size(names, 1)
    if all(names .== "")
        return String[default * string(i) for i=1:n]
    else
        for (i, n) in enumerate(names)
            names[i] = String((n == "") ? unknown : n)
        end
        return names
    end
end

function unwrap!(M::Array, dim=1)
    alldims(i) = [ n==dim ? i : (1:size(M,n)) for n in 1:ndims(M) ]
    for i = 2:size(M, dim)
        #This is a copy of slicedim from the JuliaLang but enables us to write to it
        #The code (with dim=1) is equivalent to
        # d = M[i,:,:,...,:] - M[i-1,:,...,:]
        # M[i,:,:,...,:] -= floor((d+π) / (2π)) * 2π
        d = M[alldims(i)...] - M[alldims(i-1)...]
        M[alldims(i)...] -= floor((d+π) / 2π) * 2π
    end
    return M
end

#Collect will create a copy and collect the elements
unwrap(m::AbstractArray, args...) = unwrap!(collect(m), args...)
unwrap(x::Number) = x

numpoly(G::TransferFunction) = map(numpoly, G.matrix)
denpoly(G::TransferFunction) = map(denpoly, G.matrix)

poly2vec(p::Poly) = p.a[1:end]

"""
outs = index2range(ind1, ind2)
Helper function to convert indexes with scalars to ranges. Used to avoid dropping dimensions
"""
index2range(ind1, ind2) = (index2range(ind1), index2range(ind2))
index2range{T<:Number}(ind::T) = ind:ind
index2range{T<:AbstractArray}(ind::T) = ind
index2range(ind::Colon) = ind


text:
# Get the dimensions of the problem
function ocp_dim_wrapper( Nx::Ptr{typeInt}, Nu::Ptr{typeInt}, Np::Ptr{typeInt}, Ng::Ptr{typeInt},
                          Nh::Ptr{typeInt}, NgT::Ptr{typeInt}, NhT::Ptr{typeInt}, userparam::Ptr{typeUSERPARAM} )
    if userparam == C_NULL
        error( "GRAMPC model not provided" )
    end

    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Return the problem dimension to the C code
    unsafe_store!( Nx,  model.Nx )
    unsafe_store!( Nu,  model.Nu )
    unsafe_store!( Np,  model.Np )
    unsafe_store!( Ng,  model.Ng )
    unsafe_store!( Nh,  model.Nh )
    unsafe_store!( NgT, model.NgT )
    unsafe_store!( NhT, model.NhT )
    return
end

#########################################################################################
# System dynamics
#########################################################################################

# System function f(t,x,u,p,userparam)
function ffct_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                       userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.ffct_cb, Nothing ) && return

    t  = Float64( t_in )
    x  = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    u  = unsafe_wrap( Array, u_ptr,   Int( model.Nu ) )
    p  = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    dx = unsafe_wrap( Array, out_ptr, Int( model.Nx ) )

    # Call the callback
    model.callbacks.ffct_cb( dx, x, u, p, t, model.param )
    return
end

# Jacobian df/dx multiplied by vector vec, i.e. (df/dx)^T*vec or vec^T*(df/dx)
function dfdx_vec_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, vec_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum},
                           p_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dfdx_vec_cb, Nothing ) && return

    t   = Float64( t_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    u   = unsafe_wrap( Array, u_ptr,   Int( model.Nu ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.Nx ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Nx ) )

    # Call the callback
    model.callbacks.dfdx_vec_cb( res, x, u, p, t, vec, model.param )
    return
end

# Jacobian df/du multiplied by vector vec, i.e. (df/du)^T*vec or vec^T*(df/du)
function dfdu_vec_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, vec_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum},
                           p_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dfdu_vec_cb, Nothing ) && return

    t   = Float64( t_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    u   = unsafe_wrap( Array, u_ptr,   Int( model.Nu ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.Nx ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Nu ) )

    # Call the callback
    model.callbacks.dfdu_vec_cb( res, x, u, p, t, vec, model.param )
    return
end

# Jacobian df/dp multiplied by vector vec, i.e. (df/dp)^T*vec or vec^T*(df/dp)
function dfdp_vec_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, vec_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum},
                           p_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dfdp_vec_cb, Nothing ) && return

    t   = Float64( t_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    u   = unsafe_wrap( Array, u_ptr,   Int( model.Nu ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.Nx ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Np ) )

    # Call the callback
    model.callbacks.dfdp_vec_cb( res, x, u, p, t, vec, model.param )
    return
end

#########################################################################################
# Integral cost term
#########################################################################################

# Integral cost l(t,x(t),u(t),p,xdes,udes,userparam)
function lfct_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                       xdes_ptr::Ptr{ctypeRNum}, udes_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.lfct_cb, Nothing ) && return

    t    = Float64( t_in )
    x    = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    u    = unsafe_wrap( Array, u_ptr,   Int( model.Nu ) )
    p    = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    xdes = unsafe_wrap( Array, xdes_ptr, Int( model.Nx ) )
    udes = unsafe_wrap( Array, udes_ptr, Int( model.Nu ) )

    # Call the callback
    cost = model.callbacks.lfct_cb( t, x, u, p, xdes, udes, model.param )

    unsafe_store!( out_ptr, cost )
    return
end

# Gradient dl/dx
function dldx_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                       xdes_ptr::Ptr{ctypeRNum}, udes_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dldx_cb, Nothing ) && return

    t    = Float64( t_in )
    x    = unsafe_wrap( Array, x_ptr,    Int( model.Nx ) )
    u    = unsafe_wrap( Array, u_ptr,    Int( model.Nu ) )
    p    = unsafe_wrap( Array, p_ptr,    Int( model.Np ) )
    xdes = unsafe_wrap( Array, xdes_ptr, Int( model.Nx ) )
    udes = unsafe_wrap( Array, udes_ptr, Int( model.Nu ) )
    res  = unsafe_wrap( Array, out_ptr,  Int( model.Nx ) )

    # Call the callback
    model.callbacks.dldx_cb( res, t, x, u, p, xdes, udes, model.param )
    return
end

# Gradient dl/du
function dldu_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                       xdes_ptr::Ptr{ctypeRNum}, udes_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dldu_cb, Nothing ) && return

    t    = Float64( t_in )
    x    = unsafe_wrap( Array, x_ptr,    Int( model.Nx ) )
    u    = unsafe_wrap( Array, u_ptr,    Int( model.Nu ) )
    p    = unsafe_wrap( Array, p_ptr,    Int( model.Np ) )
    xdes = unsafe_wrap( Array, xdes_ptr, Int( model.Nx ) )
    udes = unsafe_wrap( Array, udes_ptr, Int( model.Nu ) )
    res  = unsafe_wrap( Array, out_ptr,  Int( model.Nu ) )

    # Call the callback
    model.callbacks.dldu_cb( res, t, x, u, p, xdes, udes, model.param )
    return
end

# Gradient dl/dp
function dldp_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                       xdes_ptr::Ptr{ctypeRNum}, udes_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dldp_cb, Nothing ) && return

    t    = Float64( t_in )
    x    = unsafe_wrap( Array, x_ptr,    Int( model.Nx ) )
    u    = unsafe_wrap( Array, u_ptr,    Int( model.Nu ) )
    p    = unsafe_wrap( Array, p_ptr,    Int( model.Np ) )
    xdes = unsafe_wrap( Array, xdes_ptr, Int( model.Nx ) )
    udes = unsafe_wrap( Array, udes_ptr, Int( model.Nu ) )
    res  = unsafe_wrap( Array, out_ptr,  Int( model.Np ) )

    # Call the callback
    model.callbacks.dldp_cb( res, t, x, u, p, xdes, udes, model.param )
    return
end

#########################################################################################
# Terminal cost term
#########################################################################################

# Terminal cost V(T,x(T),p,xdes,userparam)
function Vfct_wrapper( out_ptr::Ptr{typeRNum}, T_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum}, xdes_ptr::Ptr{ctypeRNum},
                       userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.Vfct_cb, Nothing ) && return

    T    = Float64( T_in )
    x    = unsafe_wrap( Array, x_ptr,    Int( model.Nx ) )
    p    = unsafe_wrap( Array, p_ptr,    Int( model.Np ) )
    xdes = unsafe_wrap( Array, xdes_ptr, Int( model.Nx ) )

    # Call the callback
    cost = model.callbacks.Vfct_cb( T, x, p, xdes, model.param )

    unsafe_store!( out_ptr, cost )
    return
end

# Gradient dV/dx
function dVdx_wrapper( out_ptr::Ptr{typeRNum}, T_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum}, xdes_ptr::Ptr{ctypeRNum},
                       userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dVdx_cb, Nothing ) && return

    T    = Float64( T_in )
    x    = unsafe_wrap( Array, x_ptr,    Int( model.Nx ) )
    p    = unsafe_wrap( Array, p_ptr,    Int( model.Np ) )
    xdes = unsafe_wrap( Array, xdes_ptr, Int( model.Nx ) )
    res  = unsafe_wrap( Array, out_ptr,  Int( model.Nx ) )

    # Call the callback
    model.callbacks.dVdx_cb( res, T, x, p, xdes, model.param )
    return
end

# Gradient dV/dp
function dVdp_wrapper( out_ptr::Ptr{typeRNum}, T_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum}, xdes_ptr::Ptr{ctypeRNum},
                       userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dVdp_cb, Nothing ) && return

    T    = Float64( T_in )
    x    = unsafe_wrap( Array, x_ptr,    Int( model.Nx ) )
    p    = unsafe_wrap( Array, p_ptr,    Int( model.Np ) )
    xdes = unsafe_wrap( Array, xdes_ptr, Int( model.Nx ) )
    res  = unsafe_wrap( Array, out_ptr,  Int( model.Np ) )

    # Call the callback
    model.callbacks.dVdp_cb( res, T, x, p, xdes, model.param )
    return
end

# Gradient dV/dT
function dVdT_wrapper( out_ptr::Ptr{typeRNum}, T_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum}, xdes_ptr::Ptr{ctypeRNum},
                       userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dVdT_cb, Nothing ) && return

    T    = Float64( T_in )
    x    = unsafe_wrap( Array, x_ptr,    Int( model.Nx ) )
    p    = unsafe_wrap( Array, p_ptr,    Int( model.Np ) )
    xdes = unsafe_wrap( Array, xdes_ptr, Int( model.Nx ) )

    # Call the callback
    res = model.callbacks.dVdT_cb( T, x, p, xdes, model.param )

    unsafe_store!( out_ptr, res)
    return
end

#########################################################################################
# Equality constraints
#########################################################################################

# Equality constraints g(t,x(t),u(t),p,userparam) = 0
function gfct_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                       userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.gfct_cb, Nothing ) && return

    t   = Float64( t_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    u   = unsafe_wrap( Array, u_ptr,   Int( model.Nu ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Ng ) )

    # Call the callback
    model.callbacks.gfct_cb( res, t, x, u, p, model.param )
    return
end

# Jacobian dg/dx multiplied by vector vec, i.e. (dg/dx)^T*vec or vec^T*(dg/dx)
function dgdx_vec_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                           vec_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dgdx_cb, Nothing ) && return

    t   = Float64( t_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    u   = unsafe_wrap( Array, u_ptr,   Int( model.Nu ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.Ng ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Nx ) )

    # Call the callback
    model.callbacks.dgdx_cb( res, t, x, u, p, vec, model.param )
    return
end

# Jacobian dg/du multiplied by vector vec, i.e. (dg/du)^T*vec or vec^T*(dg/du)
function dgdu_vec_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                           vec_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dgdu_cb, Nothing ) && return

    t   = Float64( t_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    u   = unsafe_wrap( Array, u_ptr,   Int( model.Nu ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.Ng ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Nu ) )

    # Call the callback
    model.callbacks.dgdu_cb( res, t, x, u, p, vec, model.param )
    return
end

# Jacobian dg/dp multiplied by vector vec, i.e. (dg/dp)^T*vec or vec^T*(dg/dp)
function dgdp_vec_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                           vec_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dgdp_cb, Nothing ) && return

    t   = Float64( t_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    u   = unsafe_wrap( Array, u_ptr,   Int( model.Nu ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.Ng ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Np ) )

    # Call the callback
    model.callbacks.dgdp_cb( res, t, x, u, p, vec, model.param )
    return
end

#########################################################################################
# Inequality constraints
#########################################################################################

# Inequality constraints h(t,x(t),u(t),p,userparam) <= 0
function hfct_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                       userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.hfct_cb, Nothing ) && return

    t   = Float64( t_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    u   = unsafe_wrap( Array, u_ptr,   Int( model.Nu ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Nh ) )

    # Call the callback
    model.callbacks.hfct_cb( res, t, x, u, p, model.param )
    return
end

# Jacobian dh/dx multiplied by vector vec, i.e. (dh/dx)^T*vec or vec^T*(dg/dx)
function dhdx_vec_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                           vec_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dhdx_vec_cb, Nothing ) && return

    t   = Float64( t_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    u   = unsafe_wrap( Array, u_ptr,   Int( model.Nu ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.Nh ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Nx ) )

    # Call the callback
    model.callbacks.dhdx_vec_cb( res, t, x, u, p, vec, model.param )
    return
end

# Jacobian dh/du multiplied by vector vec, i.e. (dh/du)^T*vec or vec^T*(dg/du)
function dhdu_vec_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                           vec_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dhdu_vec_cb, Nothing ) && return

    t   = Float64( t_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    u   = unsafe_wrap( Array, u_ptr,   Int( model.Nu ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.Nh ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Nu ) )

    # Call the callback
    model.callbacks.dhdu_vec_cb( res, t, x, u, p, vec, model.param )
    return
end

# Jacobian dh/dp multiplied by vector vec, i.e. (dh/dp)^T*vec or vec^T*(dg/dp)
function dhdp_vec_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                           vec_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dhdp_vec_cb, Nothing ) && return

    t   = Float64( t_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    u   = unsafe_wrap( Array, u_ptr,   Int( model.Nu ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.Nh ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Np ) )

    # Call the callback
    model.callbacks.dhdp_vec_cb( res, t, x, u, p, vec, model.param )
end

#########################################################################################
# Terminal equality constraints
#########################################################################################

# Terminal equality constraints gT(T,x(T),p,userparam) = 0
function gTfct_wrapper( out_ptr::Ptr{typeRNum}, T_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.gTfct_cb, Nothing ) && return

    T   = Float64( T_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.NgT ) )

    # Call the callback
    model.callbacks.gTfct_cb( res, T, x, p, model.param )
    return
end

# Jacobian dgT/dx multiplied by vector vec, i.e. (dgT/dx)^T*vec or vec^T*(dgT/dx)
function dgTdx_vec_wrapper( out_ptr::Ptr{typeRNum}, T_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum}, vec_ptr::Ptr{ctypeRNum},
                            userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dgTdx_cb, Nothing ) && return

    T   = Float64( T_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.NgT ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Nx ) )

    # Call the callback
    model.callbacks.dgTdx_cb( res, T, x, p, vec, model.param )
    return
end

# Jacobian dgT/dp multiplied by vector vec, i.e. (dgT/dp)^T*vec or vec^T*(dgT/dp)
function dgTdp_vec_wrapper( out_ptr::Ptr{typeRNum}, T_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum}, vec_ptr::Ptr{ctypeRNum},
                            userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dgTdp_cb, Nothing ) && return

    T   = Float64( T_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.NgT ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Np ) )

    # Call the callback
    model.callbacks.dgTdp_cb( res, T, x, p, vec, model.param )
    return
end

# Jacobian dgT/dT multiplied by vector vec, i.e. (dgT/dT)^T*vec or vec^T*(dgT/dT)
function dgTdT_vec_wrapper( out_ptr::Ptr{typeRNum}, T_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum}, vec_ptr::Ptr{ctypeRNum},
                            userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dgTdT_cb, Nothing ) && return

    T   = Float64( T_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.NgT ) )

    # Call the callback
    res = model.callbacks.dgTdT_cb( T, x, p, vec, model.param )

    unsafe_store!( out_ptr, res )
    return
end

#########################################################################################
# Terminal inequality constraints
#########################################################################################

# Terminal inequality constraints hT(T,x(T),p,userparam) <= 0
function hTfct_wrapper( out_ptr::Ptr{typeRNum}, T_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.hTfct_cb, Nothing ) && return

    T   = Float64( T_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.NhT ) )

    # Call the callback
    model.callbacks.hTfct_cb( res, T, x, p, model.param )
    return
end

# Jacobian dhT/dx multiplied by vector vec, i.e. (dhT/dx)^T*vec or vec^T*(dhT/dx)
function dhTdx_vec_wrapper( out_ptr::Ptr{typeRNum}, T_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum}, vec_ptr::Ptr{ctypeRNum},
                            userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dhTdx_cb, Nothing ) && return

    T   = Float64( T_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.NhT ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Nx ) )

    # Call the callback
    model.callbacks.dhTdx_cb( res, T, x, p, vec, model.param )
    return
end

# Jacobian dhT/dp multiplied by vector vec, i.e. (dhT/dp)^T*vec or vec^T*(dhT/dp)
function dhTdp_vec_wrapper( out_ptr::Ptr{typeRNum}, T_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum}, vec_ptr::Ptr{ctypeRNum},
                            userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dhTdp_cb, Nothing ) && return

    T   = Float64( T_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.NhT ) )
    res = unsafe_wrap( Array, out_ptr, Int( model.Np ) )

    # Call the callback
    model.callbacks.dhTdp_cb( res, T, x, p, vec, model.param )
    return
end

# Jacobian dhT/dT multiplied by vector vec, i.e. (dhT/dT)^T*vec or vec^T*(dhT/dT)
function dhTdT_vec_wrapper( out_ptr::Ptr{typeRNum}, T_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum}, vec_ptr::Ptr{ctypeRNum},
                            userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    # Don't do anything if no callback is defined
    isa( model.callbacks.dhTdT_cb, Nothing ) && return

    T   = Float64( T_in )
    x   = unsafe_wrap( Array, x_ptr,   Int( model.Nx ) )
    p   = unsafe_wrap( Array, p_ptr,   Int( model.Np ) )
    vec = unsafe_wrap( Array, vec_ptr, Int( model.NhT ) )

    # Call the callback
    res = model.callbacks.dhTdT_cb( T, x, p, vec, model.param )

    unsafe_store!( out_ptr, res )
    return
end

#########################################################################################
# Functions required for semi-implicit systems
# M*dx/dt(t) = f(t0+t,x(t),u(t),p) using the solver RODAS,
#########################################################################################

# Jacobian df/dx in vector form (column-wise)
function dfdx_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                       userparam::Ptr{typeUSERPARAM} )
end

# Transpose of the Jacobian df/dx in vector form (column-wise)
function dfdxtrans_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                            userparam::Ptr{typeUSERPARAM} )
end

# Jacobian df/dt
function dfdt_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, p_ptr::Ptr{ctypeRNum},
                       userparam::Ptr{typeUSERPARAM} )
end

# Jacobian d(dH/dx)/dt
function dHdxdt_wrapper( out_ptr::Ptr{typeRNum}, t_in::ctypeRNum, x_ptr::Ptr{ctypeRNum}, u_ptr::Ptr{ctypeRNum}, vec_ptr::Ptr{ctypeRNum},
                         p_ptr::Ptr{ctypeRNum}, userparam::Ptr{typeUSERPARAM} )
end

# Mass matrix in vector form (column-wise, either banded or full matrix)
function Mfct_wrapper( out_ptr::Ptr{typeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    numElem = model.Nx * model.Nx

    # Return the columns of the mass matrix concatenated into a single vector
    out = unsafe_wrap( Array, out_ptr, Int( numElem ) )
    out = reshape( model.massMatrix, (:) )

    return
end

# Transposed mass matrix in vector form (column-wise, either banded or full matrix)
function Mtrans_wrapper( out_ptr::Ptr{typeRNum}, userparam::Ptr{typeUSERPARAM} )
    # Extract the GRAMPC problem
    model = unsafe_pointer_to_objref( userparam )::GRAMPC.Model

    numElem = model.Nx * model.Nx

    # Transpose the mass matrix then concatenate the columns into a single vector and return them
    out = unsafe_wrap( Array, out_ptr, Int( numElem ) )
    out = collect( reshape( model.massMatrix', (:) ) )

    return
end

# Julia wrapper for the problem function structure
mutable struct Cgrampcjl_callbacks
    ocp_dim::Ptr{Cvoid}

    ffct::Ptr{Cvoid}
    dfdx_vec::Ptr{Cvoid}
    dfdu_vec::Ptr{Cvoid}
    dfdp_vec::Ptr{Cvoid}

    lfct::Ptr{Cvoid}
    dldx::Ptr{Cvoid}
    dldu::Ptr{Cvoid}
    dldp::Ptr{Cvoid}

    Vfct::Ptr{Cvoid}
    dVdx::Ptr{Cvoid}
    dVdp::Ptr{Cvoid}
    dVdT::Ptr{Cvoid}

    gfct::Ptr{Cvoid}
    dgdx_vec::Ptr{Cvoid}
    dgdu_vec::Ptr{Cvoid}
    dgdp_vec::Ptr{Cvoid}

    hfct::Ptr{Cvoid}
    dhdx_vec::Ptr{Cvoid}
    dhdu_vec::Ptr{Cvoid}
    dhdp_vec::Ptr{Cvoid}

    gTfct::Ptr{Cvoid}
    dgTdx_vec::Ptr{Cvoid}
    dgTdp_vec::Ptr{Cvoid}
    dgTdT_vec::Ptr{Cvoid}

    hTfct::Ptr{Cvoid}
    dhTdx_vec::Ptr{Cvoid}
    dhTdp_vec::Ptr{Cvoid}
    dhTdT_vec::Ptr{Cvoid}

    dfdx::Ptr{Cvoid}
    dfdxtrans::Ptr{Cvoid}
    dfdt::Ptr{Cvoid}
    dHdxdt::Ptr{Cvoid}
    Mfct::Ptr{Cvoid}
    Mtrans::Ptr{Cvoid}

    function Cgrampcjl_callbacks()
        return new( @cfunction( ocp_dim_wrapper, Cvoid,
                                ( Ptr{typeInt}, Ptr{typeInt}, Ptr{typeInt}, Ptr{typeInt}, Ptr{typeInt}, Ptr{typeInt}, Ptr{typeInt}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( ffct_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM}) ),
                    @cfunction( dfdx_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dfdu_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dfdp_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( lfct_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dldx_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dldu_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dldp_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( Vfct_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dVdx_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dVdp_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dVdT_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( gfct_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dgdx_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dgdu_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dgdp_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( hfct_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dhdx_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dhdu_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dhdp_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( gTfct_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dgTdx_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dgTdp_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dgTdT_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( hTfct_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dhTdx_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dhTdp_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dhTdT_vec_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dfdx_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dfdxtrans_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dfdt_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( dHdxdt_wrapper, Cvoid,
                                ( Ptr{typeRNum}, ctypeRNum, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{ctypeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( Mfct_wrapper, Cvoid,
                                ( Ptr{typeRNum}, Ptr{typeUSERPARAM} ) ),
                    @cfunction( Mtrans_wrapper, Cvoid,
                                ( Ptr{typeRNum}, Ptr{typeUSERPARAM} ) )
            )
    end
end


text:
function result_dataframe(variable::JuMP.Containers.DenseAxisArray)

    result = Array{Float64,length(variable.axes)}(undef, length(variable.axes[2]), length(variable.axes[1]))
    # TODO: Remove this line once PowerSystems moves to Symbols
    names = Array{Symbol,1}(undef,length(variable.axes[1]))

    for t in variable.axes[2], (ix,name) in enumerate(variable.axes[1])

        result[t,ix] = JuMP.value(variable[name,t])

        names[ix] = Symbol(name)

    end

    return DataFrames.DataFrame(result, names)

end

function get_model_result(ps_m::PSI.CanonicalModel; kwargs...)

    results_df = Dict{Symbol, DataFrames.DataFrame}()

    for (k,v) in ps_m.variables

        results_df[k] = result_dataframe(v)

    end

    return results_df

end

function optimizer_log(ps_m::PSI.CanonicalModel; kwargs...)

    optimizer_log = Dict{Symbol,Any}()

    # TODO: Not all solvers can access the SolveTime() property
    #optimizer_log[:solve_time] = MOI.get(ps_m.JuMPmodel, MOI.SolveTime())
    optimizer_log[:termination_status] = JuMP.termination_status(ps_m.JuMPmodel)
    optimizer_log[:primal_status] = JuMP.primal_status(ps_m.JuMPmodel)
    optimizer_log[:dual_status] = JuMP.dual_status(ps_m.JuMPmodel)

    return optimizer_log

end


text:
# needs special attention to Cholesky, as the syntax & behavior changes in 0.4-pre

if VERSION >= v"0.5.0-dev+907"
    typealias CholType{T,S<:AbstractMatrix} Cholesky{T,S}
    chol_lower(a::Matrix) = ctranspose(chol(a))

    typealias CholTypeSparse{T} SparseArrays.CHOLMOD.Factor{T}

    chol_lower(cf::CholTypeSparse) = cf[:L]
elseif VERSION >= v"0.4.0-dev+4370"
    # error("Choleksy changes in 0.4.0-dev+4370 (PR #10862), we are still working to make it work with these changes.")

    typealias CholType{T,S<:AbstractMatrix} Cholesky{T,S}
    chol_lower(a::Matrix) = chol(a, Val{:L})

    typealias CholTypeSparse{T} SparseArrays.CHOLMOD.Factor{T}

    chol_lower(cf::CholTypeSparse) = cf[:L]
else
    typealias CholType{T,S} Cholesky{T}
    chol_lower(a::Matrix) = chol(a, :L)

    typealias CholTypeSparse{T} Base.LinAlg.CHOLMOD.CholmodFactor{T}

    chol_lower(cf::CholTypeSparse) = cf
end


text:
# Explicit
pred = [-2, -1, 0, 1, 2]
y = [0, 0, 0, 0, 0]

loss = (y - pred) .^ 2
grad = 2*pred

# AutoDiff
using Flux
using Flux.Tracker

pred = param([-2, -1, 0, 1, 2])
y = [0, 0, 0, 0, 0]
loss = (y - pred) .^ 2

# backprop
Tracker.back!(sum(loss))

pred.data
pred.grad


text:
# ASYBC.M
# This function does bias-correction based on the formula by Pope (1990)
# Luca Brugnolini 2014
# Based on Kilian

function  asybc(A,SIGMA,t,K,p)

# Input: Slope parameter estimates in compaion form,
#        variance of estimated residuals in companion form,
#        sample length, number of variables, lag order
# Output: Bias-corrected slope parameter estimates
T = t-p;
# vecSIGMAY=inv(eye((K*p)^2)-kron(A,A))*vec(SIGMA);  % Lutkepohl p.29 (2.1.39)
# SIGMAs = sparse(SIGMA[:]); As = sparse(A); eyes = sparse(eye((K*p)^2));
# vecSIGMAY = (eyes - kron(As, As))\SIGMAs; #***

vecSIGMAY = (eye((K*p)^2)-kron(A,A))\vec(SIGMA);
SIGMAY = reshape(vecSIGMAY, K*p, K*p);
I = eye(K*p, K*p);
B = A';

# There are K*p eigenvalues by construction
peigen = eigvals(A);
sumeig = zeros(K*p, K*p);

    for h = 1:K*p
        sumeig = sumeig + (peigen[h].\(I - peigen[h]*B)); #***

    end;

bias = SIGMA*(inv(I - B) + B/(I-B^2) + sumeig)/(SIGMAY); #***

Abias = - bias/T;

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

bcstab = 9;         # Arbitrary default value > 1
delta = 1;            # Adjustment factor

while bcstab >= 1

# Adjust bias-correction proportionately
global bcA = A-delta*Abias;
bcmod = abs(eigvals(bcA));

if  any(bcmod.>= 1)
    bcstab = 1;

else
    bcstab = 0;

end

delta += - 0.01;

if delta <= 0
   bcstab = 0;

end

end

return bcA

end


text:
mutable struct SumCache{Ta,Tb,Tc}
    σᵢₙ::Ta
    out::Tc
    δℒ::Tb
    valid::Bool
end

cache(::typeof(sum_autobatch), arr_T, in_T, in_sz) =
    SumCache(similar(arr_T, in_T, in_sz...),
              zero(in_T),
              similar(arr_T, in_T, prod(in_sz)),
              false)

layer_out_type_size(::typeof(sum_autobatch), in_T ,in_sz) = in_T, (1,)

function sum_autobatch(c::SumCache, x::AbstractVector)
    σ = copyto!(c.σᵢₙ, x)

    c.out = sum(x)
    return c.out
end

function backprop(∇, l::typeof(sum_autobatch), c::SumCache, δℒ::Number)
    c.δℒ .= δℒ

    return c.δℒ
end


text:
using Koma, JLD2
#Scanner
sys = Scanner()
sys.ADC_Δt = 4e-6
#Seq
α, T = π/2, 1e-3
B1 = α / (2π*γ*T) # α =  2π γ B1 T -> B1 = α / (2π γ T)
EX = PulseDesigner.RF_hard(B1,T,sys)
FOV, N = 25.6e-2, 101
EPI,_,_,_ = PulseDesigner.EPI_base(FOV,N,sys)
TE = 50e-3
d1 = TE - dur(EPI)/2 - dur(EX)
D1 = Delay(d1) #You can obtain the duration of the EPI with dur(EPI)
R = rotz(0.) #Play with the rotation of the acquisition!
# For a GE acquisition
seq = EX + D1 + R*EPI
#Save seq
println(seq)
@save "./myEPI_GE.seq" seq
# For a SE, the delays must refocuse at k=0 (half the EPI)
REF = 2.0im*EX #twice the B1, 90 deg -> 180 deg
halfTE = TE / 2
d1 = halfTE - dur(EX)
d2 = halfTE - dur(EPI)/2 - dur(REF)
D1 = Delay(d1)
D2 = Delay(d2)
seq = EX + D1 + REF + D2 + R*EPI
#Save seq
println(seq)
@save "./myEPI_SE.seq" seq
#Go to SpinLab GUI, loag "myEPI_GE.seq", press "Run"