text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.cluster;
#
# import typedb.client.common.exception.TypeDBClientException;
# import typedb.protocol.ClusterDatabaseProto;
# import io.grpc.StatusRuntimeException;
# import org.slf4j.Logger;
# import org.slf4j.LoggerFactory;
#
# import java.util.ArrayList;
# import java.util.List;
#
# import static typedb.client.common.exception.ErrorMessage.Client.CLUSTER_REPLICA_NOT_PRIMARY;
# import static typedb.client.common.exception.ErrorMessage.Client.CLUSTER_UNABLE_TO_CONNECT;
# import static typedb.client.common.exception.ErrorMessage.Client.UNABLE_TO_CONNECT;
# import static typedb.client.common.exception.ErrorMessage.Internal.UNEXPECTED_INTERRUPTION;
# import static typedb.client.common.rpc.RequestBuilder.Cluster.DatabaseManager.getReq;
#
# abstract class FailsafeTask<RESULT> {
#
#     private static final Logger LOG = LoggerFactory.getLogger(FailsafeTask.class);
#     private static final int PRIMARY_REPLICA_TASK_MAX_RETRIES = 10;
#     private static final int FETCH_REPLICAS_MAX_RETRIES = 10;
#     private static final int WAIT_FOR_PRIMARY_REPLICA_SELECTION_MS = 2000;
#
#     private final ClusterClient client;
#     private final String database;
#
#     FailsafeTask(ClusterClient client, String database) {
#         this.client = client;
#         this.database = database;
#     }
#
#     abstract RESULT run(ClusterDatabase.Replica replica);
#
#     RESULT rerun(ClusterDatabase.Replica replica) {
#         return run(replica);
#     }
#
#     RESULT runPrimaryReplica() {
#         if (!client.databaseByName().containsKey(database)
#                 || !client.databaseByName().get(database).primaryReplica().isPresent()) {
#             seekPrimaryReplica();
#         }
#         ClusterDatabase.Replica replica = client.databaseByName().get(database).primaryReplica().get();
#         int retries = 0;
#         while (true) {
#             try {
#                 return retries == 0 ? run(replica) : rerun(replica);
#             } catch (TypeDBClientException e) {
#                 if (CLUSTER_REPLICA_NOT_PRIMARY.equals(e.getErrorMessage())
#                         || UNABLE_TO_CONNECT.equals(e.getErrorMessage())) {
#                     LOG.debug("Unable to open a session or transaction, retrying in 2s...", e);
#                     waitForPrimaryReplicaSelection();
#                     replica = seekPrimaryReplica();
#                 } else throw e;
#             }
#             if (++retries > PRIMARY_REPLICA_TASK_MAX_RETRIES) throw clusterNotAvailableException();
#         }
#     }
#
#     RESULT runAnyReplica() {
#         ClusterDatabase databaseClusterRPC = client.databaseByName().get(database);
#         if (databaseClusterRPC == null) databaseClusterRPC = fetchDatabaseReplicas();
#
#         // Try the preferred secondary replica first, then go through the others
#         List<ClusterDatabase.Replica> replicas = new ArrayList<>();
#         replicas.add(databaseClusterRPC.preferredReplica());
#         for (ClusterDatabase.Replica replica : databaseClusterRPC.replicas()) {
#             if (!replica.isPreferred()) replicas.add(replica);
#         }
#
#         int retries = 0;
#         for (ClusterDatabase.Replica replica : replicas) {
#             try {
#                 return retries == 0 ? run(replica) : rerun(replica);
#             } catch (TypeDBClientException e) {
#                 if (UNABLE_TO_CONNECT.equals(e.getErrorMessage())) {
#                     LOG.debug("Unable to open a session or transaction to " + replica.id() +
#                                       ". Attempting next replica.", e);
#                 } else {
#                     throw e;
#                 }
#             }
#             retries++;
#         }
#         throw clusterNotAvailableException();
#     }
#
#     private ClusterDatabase.Replica seekPrimaryReplica() {
#         int retries = 0;
#         while (retries < FETCH_REPLICAS_MAX_RETRIES) {
#             ClusterDatabase databaseClusterRPC = fetchDatabaseReplicas();
#             if (databaseClusterRPC.primaryReplica().isPresent()) {
#                 return databaseClusterRPC.primaryReplica().get();
#             } else {
#                 waitForPrimaryReplicaSelection();
#                 retries++;
#             }
#         }
#         throw clusterNotAvailableException();
#     }
#
#     private ClusterDatabase fetchDatabaseReplicas() {
#         for (String serverAddress : client.clusterMembers()) {
#             try {
#                 LOG.debug("Fetching replica info from {}", serverAddress);
#                 ClusterDatabaseProto.ClusterDatabaseManager.Get.Res res = client.stub(serverAddress).databasesGet(getReq(database));
#                 ClusterDatabase databaseClusterRPC = ClusterDatabase.of(res.getDatabase(), client.databases());
#                 client.databaseByName().put(database, databaseClusterRPC);
#                 return databaseClusterRPC;
#             } catch (TypeDBClientException e) {
#                 if (e.getErrorMessage().equals(UNABLE_TO_CONNECT)) {
#                     LOG.debug("Failed to fetch replica info for database '" + database + "' from " +
#                             serverAddress + ". Attempting next server.", e);
#                 } else {
#                     throw e;
#                 }
#             }
#         }
#         throw clusterNotAvailableException();
#     }
#
#     private void waitForPrimaryReplicaSelection() {
#         try {
#             Thread.sleep(WAIT_FOR_PRIMARY_REPLICA_SELECTION_MS);
#         } catch (InterruptedException e) {
#             throw new TypeDBClientException(UNEXPECTED_INTERRUPTION);
#         }
#     }
#
#     private TypeDBClientException clusterNotAvailableException() {
#         return new TypeDBClientException(CLUSTER_UNABLE_TO_CONNECT, String.join(",", client.clusterMembers()));
#     }
# }


text:
#=
# FV solver
=#
module FVSolver

export solve, nonlinear_solve

import ..Femshop: JULIA, CPP, MATLAB, DENDRO, HOMG, CUSTOM_GEN_TARGET,
                SQUARE, IRREGULAR, UNIFORM_GRID, TREE, UNSTRUCTURED,
                CG, DG, HDG, FV,
                NODAL, MODAL, CELL, LEGENDRE, UNIFORM, GAUSS, LOBATTO,
                NONLINEAR_NEWTON, NONLINEAR_SOMETHING,
                EULER_EXPLICIT, EULER_IMPLICIT, CRANK_NICHOLSON, RK4, LSRK4, ABM4,
                DEFAULT_SOLVER, PETSC,
                VTK, RAW_OUTPUT, CUSTOM_OUTPUT,
                DIRICHLET, NEUMANN, ROBIN, NO_BC, FLUX,
                MSH_V2, MSH_V4,
                SCALAR, VECTOR, TENSOR, SYM_TENSOR, VAR_ARRAY,
                LHS, RHS,
                LINEMESH, QUADMESH, HEXMESH
import ..Femshop: log_entry, printerr
import ..Femshop: config, prob, variables, mesh_data, grid_data, refel, time_stepper, elemental_order, genfunctions, indexers
import ..Femshop: Variable, Coefficient, GenFunction
import ..Femshop: GeometricFactors, geo_factors, geometric_factors, geometric_factors_face, build_deriv_matrix
# specific to FV
import ..Femshop: FVInfo, fv_info, FV_cell_to_node, FV_node_to_cell
import ..Femshop: ParentMaps, parent_maps, fv_grid, fv_refel, fv_geo_factors, build_local_patch

using LinearAlgebra, SparseArrays

include("fv_boundary.jl");

function linear_solve(var, source_lhs, source_rhs, flux_lhs, flux_rhs, stepper=nothing, assemble_func=nothing)
    # If more than one variable
    if typeof(var) <: Array
        # multiple variables being solved for simultaneously
        dofs_per_node = 0;
        dofs_per_loop = 0;
        for vi=1:length(var)
            dofs_per_loop += length(var[vi].symvar);
            dofs_per_node += var[vi].total_components;
        end
    else
        # one variable
        dofs_per_loop = length(var.symvar);
        dofs_per_node = var.total_components;
    end
    if fv_grid === nothing
        nel = size(grid_data.loc2glb, 2);
        nfaces = size(grid_data.face2element, 2);
    else
        nel = size(fv_grid.loc2glb, 2);
        nfaces = size(fv_grid.face2element, 2);
    end

    Nn = dofs_per_node * nel;
    Nf = dofs_per_node * nfaces

    # Allocate arrays that will be used by assemble
    # These vectors will hold the integrated values(one per cell).
    # They will later be combined.
    sourcevec = zeros(Nn);
    fluxvec = zeros(Nn);
    facefluxvec = zeros(Nf);
    face_done = zeros(Int, nfaces); # Increment when the corresponding flux value is computed.
    allocated_vecs = [sourcevec, fluxvec, facefluxvec, face_done];

    if prob.time_dependent && !(stepper === nothing)
        log_entry("Beginning "*string(stepper.Nsteps)*" time steps.");
        t = 0;
        sol = get_var_vals(var);

        # allocate storage used by steppers
        if stepper.type == LSRK4
            tmppi = zeros(size(sol));
            tmpki = zeros(size(sol));
        elseif stepper.type == RK4
            tmpki = zeros(length(sol), stepper.stages);
        end

        start_t = Base.Libc.time();
        last2update = 0;
        last10update = 0;
        print("Time stepping progress(%): 0");
        for i=1:stepper.Nsteps
            if stepper.stages > 1
                # LSRK4 is a special case, low storage
                if stepper.type == LSRK4
                    # Low storage RK4:
                    # p0 = u
                    #   ki = ai*k(i-1) + dt*f(p(i-1), t+ci*dt)
                    #   pi = p(i-1) + bi*ki
                    # u = p5

                    tmppi = get_var_vals(var, tmppi);
                    tmpki = zeros(size(sol));
                    for rki=1:stepper.stages
                        rktime = t + stepper.c[rki]*stepper.dt;
                        # p(i-1) is currently in u

                        sol = assemble(var, source_lhs, source_rhs, flux_lhs, flux_rhs, allocated_vecs, dofs_per_node, dofs_per_loop, stime, stepper.dt, assemble_loops=assemble_func);

                        if rki == 1 # because a1 == 0
                            tmpki = stepper.dt .* sol;
                        else
                            tmpki = stepper.a[rki].*tmpki + stepper.dt.*sol;
                        end
                        tmppi = tmppi + stepper.b[rki].*tmpki

                        place_sol_in_vars(var, tmppi, stepper);
                    end

                else
                    # Explicit multi-stage methods:
                    # x = x + dt*sum(bi*ki)
                    # ki = rhs(t+ci*dt, x+dt*sum(aij*kj)))   j < i

                    # will hold the final result
                    sol = get_var_vals(var, sol);
                    # will be placed in var.values for each stage
                    tmpvals = sol;
                    for stage=1:stepper.stages
                        stime = t + stepper.c[stage]*stepper.dt;

                        tmpki[:,stage] = assemble(var, source_lhs, source_rhs, flux_lhs, flux_rhs, allocated_vecs, dofs_per_node, dofs_per_loop, stime, stepper.dt, assemble_loops=assemble_func);

                        tmpvals = sol;
                        for j=1:(stage-1)
                            if stepper.a[stage, j] > 0
                                tmpvals += stepper.dt * stepper.a[stage, j] .* tmpki[:,j];
                            end
                        end

                        place_sol_in_vars(var, tmpvals, stepper);
                    end
                    for stage=1:stepper.stages
                        sol += stepper.dt * stepper.b[stage] .* tmpki[:, stage];
                    end
                    place_sol_in_vars(var, sol, stepper);
                end

            elseif stepper.type == EULER_EXPLICIT
                sol = sol .+ stepper.dt .* assemble(var, source_lhs, source_rhs, flux_lhs, flux_rhs, allocated_vecs, dofs_per_node, dofs_per_loop, t, stepper.dt, assemble_loops=assemble_func);

                place_sol_in_vars(var, sol, stepper);

            else
                printerr("Only explicit time steppers for FV are ready. TODO")
                return sol;
            end

            # ########## uncomment to return after one time step
            # return sol
            # ############

            t += stepper.dt;

            progressPercent = Int(floor(i*100.0/stepper.Nsteps));
            if progressPercent - last2update >= 2
                last2update = progressPercent;
                if progressPercent - last10update >= 10
                    print(progressPercent);
                    last10update = progressPercent;
                else
                    print(".");
                end
            end
        end
        println("");
        end_t = Base.Libc.time();

        log_entry("Stepping took "*string(end_t-start_t)*" seconds.");
        #log_entry("Stepping took "*string(end_t-start_t)*" seconds. ("*string(assemble_t)*" for assembly, "*string(linsolve_t)*" for linear solve)");
        #display(sol);
		# outfile = "linear_sol.txt"
		# open(outfile, "w") do f
  		# 	for ii in sol
    	# 		println(f, ii)
  		# 	end
		# end # the file f is automatically closed after this block finishes
        return sol;

    else
        # Does it make any sense to do this for time-independent problems?
        return sol;
    end
end

# macro loops(indices, ranges, content)
#     n = length(indices.args);
#     if n == 1
#         this_ind = indices.args[1];
#         this_range = ranges.args[1];
#         return esc(quote
#             for $this_ind in $this_range
#                 $content
#             end
#         end)
#     else
#         this_ind = indices.args[1];
#         this_range = ranges.args[1];
#         sub_ind = copy(indices);
#         sub_ind.args = sub_ind.args[2:end];
#         sub_range = copy(ranges);
#         sub_range.args = sub_range.args[2:end];
#         return esc(:(@loops([$this_ind], [$this_range], @loops($sub_ind, $sub_range, $content))));
#     end
# end

#

function assemble(var, source_lhs, source_rhs, flux_lhs, flux_rhs, allocated_vecs, dofs_per_node=1, dofs_per_loop=1, t=0, dt=0; assemble_loops=nothing)
    # If an assembly loop function was provided, use it
    if !(assemble_loops === nothing)
        return assemble_loops.func(var, source_lhs, source_rhs, flux_lhs, flux_rhs, allocated_vecs, dofs_per_node, dofs_per_loop, t, dt);
    end
    # If parent maps were created for high order flux, use a slightly different function
    if !(parent_maps === nothing)
        return assemble_using_parent_child(var, source_lhs, source_rhs, flux_lhs, flux_rhs, allocated_vecs, dofs_per_node, dofs_per_loop, t, dt);
    end

    nel = size(grid_data.loc2glb, 2);

    # Label things that were allocated externally
    sourcevec = allocated_vecs[1];
    fluxvec = allocated_vecs[2];
    facefluxvec = allocated_vecs[3];
    face_done = allocated_vecs[4];

    face_done .= 0;

    # Elemental loop
    for ei=1:nel
        eid = elemental_order[ei]; # The index of this element
        # Zero the result vectors for this element
        sourcevec[((eid-1)*dofs_per_node + 1):(eid*dofs_per_node)] .= 0;
        fluxvec[((eid-1)*dofs_per_node + 1):(eid*dofs_per_node)] .= 0;

        ##### Source integrated over the cell #####
        # Compute RHS volume integral
        if !(source_rhs === nothing)
            #sourceargs = prepare_args(var, eid, 0, RHS, "volume", t, dt); # (var, e, nodex, loc2glb, refel, detj, J, t, dt)
            sourceargs = (var, eid, 0, grid_data, geo_factors, fv_info, refel, t, dt);
            source = source_rhs.func(sourceargs) ./ geo_factors.volume[eid];
            # Add to global source vector
            sourcevec[((eid-1)*dofs_per_node + 1):(eid*dofs_per_node)] = source;
        end

        ##### Flux integrated over the faces #####
        # Loop over this element's faces.
        for i=1:refel.Nfaces
            fid = grid_data.element2face[i, eid];

            if !(flux_rhs === nothing)
                if face_done[fid] == 0
                    face_done[fid] = 1; # Possible race condition, but in the worst case it will be computed twice.

                    #fluxargs = prepare_args(var, eid, fid, RHS, "surface", t, dt); #(var, (e, neighbor), refel, vol_loc2glb, nodex, cellx, frefelind, facex, face2glb, normal, fdetj, face_area, (J, vol_J_neighbor), t, dt);
                    fluxargs = (var, eid, fid, grid_data, geo_factors, fv_info, refel, t, dt);
                    flux = flux_rhs.func(fluxargs) .* geo_factors.area[fid];
                    # Add to global flux vector for faces
                    facefluxvec[((fid-1)*dofs_per_node + 1):(fid*dofs_per_node)] = flux;
                    # Combine all flux for this element
                    fluxvec[((eid-1)*dofs_per_node + 1):(eid*dofs_per_node)] += flux ./ geo_factors.volume[eid];

                else
                    # This flux has either been computed or is being computed by another thread.
                    # The state will need to be known before paralellizing, but for now assume it's complete.
                    fluxvec[((eid-1)*dofs_per_node + 1):(eid*dofs_per_node)] -= facefluxvec[((fid-1)*dofs_per_node + 1):(fid*dofs_per_node)] ./ geo_factors.volume[eid];
                end
            end

            # Boundary conditions are applied to flux
            fbid = grid_data.facebid[fid]; # BID of this face
            if fbid > 0
                facex = grid_data.allnodes[:, grid_data.face2glb[:,1,fid]];  # face node coordinates

                if typeof(var) <: Array
                    dofind = 0;
                    for vi=1:length(var)
                        for compo=1:length(var[vi].symvar)
                            dofind = dofind + 1;
                            if prob.bc_type[var[vi].index, fbid] == NO_BC
                                # do nothing
                            elseif prob.bc_type[var[vi].index, fbid] == FLUX
                                # compute the value and add it to the flux directly
                                # Qvec = (refel.surf_wg[grid_data.faceRefelInd[1,fid]] .* geo_factors.face_detJ[fid])' * (refel.surf_Q[grid_data.faceRefelInd[1,fid]])[:, refel.face2local[grid_data.faceRefelInd[1,fid]]]
                                # Qvec = Qvec ./ geo_factors.area[fid];
                                # bflux = FV_flux_bc_rhs_only(prob.bc_func[var[vi].index, fbid][compo], facex, Qvec, t, dofind, dofs_per_node) .* geo_factors.area[fid];
                                bflux = FV_flux_bc_rhs_only_simple(prob.bc_func[var[vi].index, fbid][compo], facex, t) .* geo_factors.area[fid];

                                fluxvec[(eid-1)*dofs_per_node + dofind] += (bflux - facefluxvec[(fid-1)*dofs_per_node + dofind]) ./ geo_factors.volume[eid];
                                facefluxvec[(fid-1)*dofs_per_node + dofind] = bflux;
                            else
                                printerr("Unsupported boundary condition type: "*prob.bc_type[var[vi].index, fbid]);
                            end
                        end
                    end
                else
                    for d=1:dofs_per_node
                        dofind = d;
                        if prob.bc_type[var.index, fbid] == NO_BC
                            # do nothing
                        elseif prob.bc_type[var.index, fbid] == FLUX
                            # compute the value and add it to the flux directly
                            # Qvec = (refel.surf_wg[grid_data.faceRefelInd[1,fid]] .* geo_factors.face_detJ[fid])' * (refel.surf_Q[grid_data.faceRefelInd[1,fid]])[:, refel.face2local[grid_data.faceRefelInd[1,fid]]]
                            # Qvec = Qvec ./ geo_factors.area[fid];
                            # bflux = FV_flux_bc_rhs_only(prob.bc_func[var.index, fbid][d], facex, Qvec, t, dofind, dofs_per_node) .* geo_factors.area[fid];
                            bflux = FV_flux_bc_rhs_only_simple(prob.bc_func[var.index, fbid][d], facex, t) .* geo_factors.area[fid];

                            fluxvec[(eid-1)*dofs_per_node + dofind] += (bflux - facefluxvec[(fid-1)*dofs_per_node + dofind]) ./ geo_factors.volume[eid];
                            facefluxvec[(fid-1)*dofs_per_node + dofind] = bflux;
                        else
                            printerr("Unsupported boundary condition type: "*prob.bc_type[var.index, fbid]);
                        end
                    end
                end
            end# BCs

        end# face loop
    end# element loop

    return sourcevec + fluxvec;
end

function assemble_using_parent_child(var, source_lhs, source_rhs, flux_lhs, flux_rhs, allocated_vecs, dofs_per_node=1, dofs_per_loop=1, t=0, dt=0)
    debug = false;

    # Label things that were allocated externally
    sourcevec = allocated_vecs[1];
    fluxvec = allocated_vecs[2];
    facefluxvec = allocated_vecs[3];
    face_done = allocated_vecs[4];

    face_done .= 0;
    fluxvec .= 0; # need to do this here

    # Loop over parent elements
    # Since flux depends only on this and nearest neighbor parent cells, it is safe
    # to partition between parents using only the nearest parents as ghosts.
    Nparents = size(parent_maps.parent2child,2);
    nchildren = size(parent_maps.parent2child,1); # (assumes one element type)
    npfaces = size(grid_data.element2face,1); # Outer faces of parent. Parents are elements in grid_data (assumes one element type)
    ncfaces = size(parent_maps.parent2face,1); # All faces of children in a parent (assumes one element type)
    dim = size(grid_data.allnodes,1);

    for parentid=1:Nparents
        # Get the cells in the local patch
        patch_cells = parent_maps.patches[:, parentid];

        # Split the work into two parts:
        # - source terms that don't depend on neighbors and can be done independently (loop over elements)
        # - flux terms that depend on neighbors (loop over faces)

        # The loop over elements for source terms.
        for childid = 1:nchildren
            eid = parent_maps.parent2child[childid, parentid];
            # Zero the result vectors for this element
            sourcevec[((eid-1)*dofs_per_node + 1):(eid*dofs_per_node)] .= 0;

            # Use source_rhs for RHS volume integral
            if !(source_rhs === nothing)
                sourceargs = (var, eid, 0, fv_grid, fv_geo_factors, fv_info, fv_refel, t, dt);
                source = source_rhs.func(sourceargs; patch_cells=patch_cells) ./ fv_geo_factors.volume[eid];
                # Add to global source vector
                sourcevec[((eid-1)*dofs_per_node + 1):(eid*dofs_per_node)] = source;
            end
        end

        # Loop over this parent's faces for flux terms.
        for pfaceid = 1:ncfaces
            fid = parent_maps.parent2face[pfaceid, parentid];
            els = fv_grid.face2element[:,fid];
            # # The eid passed to the flux func. will be one inside the parent.
            # if els[1] in parent_maps.parent2child[:,parentid]
            #     eid = els[1];
            #     neighbor = els[2];
            # else
            #     eid = els[2];
            #     neighbor = els[1];
            # end
            eid = els[1];
            neighbor = els[2];

            if debug && parentid < 3
                println("parent: "*string(parentid)*", fid: "*string(fid)*"("*string(pfaceid)*"), eid/neigh: "*string(eid)*"/"*string(neighbor));
                println("patch: "*string(patch_cells));
            end

            if !(flux_rhs === nothing)
                if face_done[fid] == 0
                    face_done[fid] = 1; # Possible race condition, but in the worst case it will be computed twice.\

                    (left_cells, right_cells) = get_left_right_cells(patch_cells, pfaceid, parent_maps, dim);

                    # leftness and rightness depends on the normal vector direction
                    # normal points left to right.
                    leftcell = fv_grid.face2element[1,fid];
                    if leftcell == right_cells[1]
                        #oops, swap them
                        tmp = left_cells;
                        left_cells = right_cells;
                        right_cells = tmp;
                    end

                    # remove zero values
                    left_cells = remove_zero_cells(left_cells);
                    right_cells = remove_zero_cells(right_cells);

                    #fluxargs = prepare_args(var, eid, fid, RHS, "surface", t, dt); #(var, (e, neighbor), refel, vol_loc2glb, nodex, cellx, frefelind, facex, face2glb, normal, fdetj, face_area, (J, vol_J_neighbor), t, dt);
                    fluxargs = (var, eid, fid, fv_grid, fv_geo_factors, fv_info, fv_refel, t, dt);
                    flux = flux_rhs.func(fluxargs; left_cells=left_cells, right_cells=right_cells) .* fv_geo_factors.area[fid];

                    # Add to global flux vector for faces
                    facefluxvec[((fid-1)*dofs_per_node + 1):(fid*dofs_per_node)] = flux;
                    # Contribute to flux for elements on both sides
                    fluxvec[((eid-1)*dofs_per_node + 1):(eid*dofs_per_node)] += flux ./ fv_geo_factors.volume[eid];
                    if neighbor > 0 # not a boundary, so the other side exists
                        fluxvec[((neighbor-1)*dofs_per_node + 1):(neighbor*dofs_per_node)] -= flux ./ fv_geo_factors.volume[neighbor];
                    end

                    if debug && parentid < 3
                        println("cells left: "*string(left_cells)*", right: "*string(right_cells));
                        println("flux: "*string(flux));
                        println("contribute "*string(flux ./ fv_geo_factors.volume[eid])*" to "*string(((eid-1)*dofs_per_node + 1):(eid*dofs_per_node)));
                        if neighbor > 0
                            println("contribute "*string(-flux ./ fv_geo_factors.volume[neighbor])*" to "*string(((neighbor-1)*dofs_per_node + 1):(neighbor*dofs_per_node)));
                        end
                        println("fluxvec "*string(fluxvec[1:8]));
                        println("");
                    end

                else
                    # This flux has either been computed or is being computed by another thread.
                    # The state will need to be known before paralellizing, but for now assume it's complete.
                    # fluxvec[((eid-1)*dofs_per_node + 1):(eid*dofs_per_node)] .-= facefluxvec[((fid-1)*dofs_per_node + 1):(fid*dofs_per_node)] ./ fv_geo_factors.volume[eid];
                end
            end

            # Boundary conditions are applied to flux
            fbid = fv_grid.facebid[fid]; # BID of this face
            if fbid > 0
                if debug
                    println("applying bdry cond"*string(fbid)*" to "*string(fid));
                end
                facex = fv_grid.allnodes[:, fv_grid.face2glb[:,1,fid]];  # face node coordinates

                if typeof(var) <: Array
                    dofind = 0;
                    for vi=1:length(var)
                        for compo=1:length(var[vi].symvar)
                            dofind = dofind + 1;
                            if prob.bc_type[var[vi].index, fbid] == NO_BC
                                # do nothing
                            elseif prob.bc_type[var[vi].index, fbid] == FLUX
                                # compute the value and add it to the flux directly
                                Qvec = (fv_refel.surf_wg[fv_grid.faceRefelInd[1,fid]] .* fv_geo_factors.face_detJ[fid])' * (fv_refel.surf_Q[fv_grid.faceRefelInd[1,fid]])[:, fv_refel.face2local[fv_grid.faceRefelInd[1,fid]]]
                                Qvec = Qvec ./ fv_geo_factors.area[fid];
                                bflux = FV_flux_bc_rhs_only(prob.bc_func[var[vi].index, fbid][compo], facex, Qvec, t, dofind, dofs_per_node) .* fv_geo_factors.area[fid];

                                fluxvec[(eid-1)*dofs_per_node + dofind] += (bflux - facefluxvec[(fid-1)*dofs_per_node + dofind]) ./ fv_geo_factors.volume[eid];
                                facefluxvec[(fid-1)*dofs_per_node + dofind] = bflux;
                            else
                                printerr("Unsupported boundary condition type: "*prob.bc_type[var[vi].index, fbid]);
                            end
                        end
                    end
                else
                    for d=1:dofs_per_node
                        dofind = d;
                        if prob.bc_type[var.index, fbid] == NO_BC
                            # do nothing
                        elseif prob.bc_type[var.index, fbid] == FLUX
                            # compute the value and add it to the flux directly
                            Qvec = (fv_refel.surf_wg[fv_grid.faceRefelInd[1,fid]] .* fv_geo_factors.face_detJ[fid])' * (fv_refel.surf_Q[fv_grid.faceRefelInd[1,fid]])[:, fv_refel.face2local[fv_grid.faceRefelInd[1,fid]]]
                            Qvec = Qvec ./ fv_geo_factors.area[fid];
                            bflux = FV_flux_bc_rhs_only(prob.bc_func[var.index, fbid][d], facex, Qvec, t, dofind, dofs_per_node) .* fv_geo_factors.area[fid];

                            fluxvec[(eid-1)*dofs_per_node + dofind] += (bflux - facefluxvec[(fid-1)*dofs_per_node + dofind]) ./ fv_geo_factors.volume[eid];
                            facefluxvec[(fid-1)*dofs_per_node + dofind] = bflux;
                        else
                            printerr("Unsupported boundary condition type: "*prob.bc_type[var.index, fbid]);
                        end
                    end
                end
            end# BCs
        end# face loop
    end# element loop
    if debug
        println("fluxvec: "*string(fluxvec[1:10]));
    end
    return sourcevec + fluxvec;
end

# Returns a copy of a with zeros removed.
function remove_zero_cells(a)
    b = similar(a);
    nnz = 0;
    for ai in a
        if !(ai == 0)
            nnz += 1;
            b[nnz] = ai;
        end
    end
    b = b[1:nnz];

    return b;
end

# Returns arrays of left and right cells
function get_left_right_cells(patch, face, maps, dim)
    # Provide all of the cells that can be used for this face in a particular order.
    # Make two arrays, one for left one for right, starting from the nearest neighbor.
    left_cells = [];
    right_cells = [];
    if dim == 1
        left_cell_table_1d = [
            [ # children=1
                [2],
                [1, 2]
            ],
            [ # children=2
                [4, 3],
                [1, 4, 3],
                [2, 1, 4, 3]
            ],
            [ # children=3
                [6, 5, 4],
                [1, 6, 5, 4],
                [2, 1, 6, 5, 4],
                [3, 2, 1, 6, 5, 4]
            ],
            [ # children=4
                [8, 7, 6, 5],
                [1, 8, 7, 6, 5],
                [2, 1, 8, 7, 6, 5],
                [3, 2, 1, 8, 7, 6, 5],
                [4, 3, 2, 1, 8, 7, 6, 5],
            ]
        ]
        right_cell_table_1d = [
            [ # children=1
                [1, 3],
                [3]
            ],
            [ # children=2
                [1, 2, 5, 6],
                [2, 5, 6],
                [5, 6]
            ],
            [ # children=3
                [1, 2, 3, 7, 8, 9],
                [2, 3, 7, 8, 9],
                [3, 7, 8, 9],
                [7, 8, 9]
            ],
            [ # children=4
                [1, 2, 3, 4, 9, 10, 11, 12],
                [2, 3, 4, 9, 10, 11, 12],
                [3, 4, 9, 10, 11, 12],
                [4, 9, 10, 11, 12],
                [9, 10, 11, 12]
            ]
        ]
        nchildren = size(maps.parent2child,1);
        left_cells = left_cell_table_1d[nchildren][face];
        right_cells = right_cell_table_1d[nchildren][face];

    elseif dim == 2
        if size(maps.parent2neighbor,1) == 3 # triangles
            # Triangle parents have 9 faces, patches have 16 cells
            # Here Left means toward the center of the central parent
            left_cell_table_triangle = [
                [1, 4, 3, 15, 16, 13, 2, ],
                [2, 4, 3, 9, 12, 11, 1],
                [2, 4, 1, 7, 8, 5, 3],
                [3, 4, 1, 13, 16, 15, 2],
                [3, 4, 2, 11, 12, 9, 1],
                [1, 4, 2, 5, 8, 7, 3],
                [4, 2, 3, 9, 11, 12],
                [4, 1, 3, 15, 13, 16],
                [4, 1, 2, 5, 7, 8]
            ]
            right_cell_table_triangle = [
                [5, 8, 6, 7],
                [7, 8, 6, 5],
                [9, 12, 10, 11],
                [11, 12, 10, 9],
                [13, 16, 14, 15],
                [15, 16, 14, 13],
                [1, 5, 15, 8, 16],
                [2, 7, 9, 8, 12],
                [3, 11, 13, 12, 16]
            ]
            left_cells = left_cell_table_triangle[face];
            right_cells = right_cell_table_triangle[face];
        else # quads
            #TODO
        end

    elseif dim == 3
        #TODO
    end

    return (patch[left_cells], patch[right_cells]);
end

# Insert the single dof into the greater vector
function insert_linear!(b, bel, glb, dof, Ndofs)
    # group nodal dofs
    for d=1:length(dof)
        ind = glb.*Ndofs .- (Ndofs-dof[d]);
        ind2 = ((d-1)*length(glb)+1):(d*length(glb));

        b[ind] = b[ind] + bel[ind2];
    end
end

function place_sol_in_vars(var, sol, stepper)
    # place the values in the variable value arrays
    if typeof(var) <: Array
        tmp = 0;
        totalcomponents = 0;
        for vi=1:length(var)
            totalcomponents = totalcomponents + var[vi].total_components;
        end
        for vi=1:length(var)
            components = var[vi].total_components;
            for compi=1:components
                if stepper.type == EULER_EXPLICIT
                    var[vi].values[compi,:] = sol[(compi+tmp):totalcomponents:end];
                else
                    var[vi].values[compi,:] = sol[(compi+tmp):totalcomponents:end];
                end
                tmp = tmp + 1;
            end
        end
    else
        components = var.total_components;
        for compi=1:components
            if stepper.type == EULER_EXPLICIT
                var.values[compi,:] = sol[compi:components:end];
            else
                var.values[compi,:] = sol[compi:components:end];
            end
        end
    end
end

function get_var_vals(var, vect=nothing)
    # place the variable values in a vector
    if typeof(var) <: Array
        tmp = 0;
        totalcomponents = 0;
        for vi=1:length(var)
            totalcomponents = totalcomponents + var[vi].total_components;
        end
        if vect === nothing
            vect = zeros(totalcomponents * length(var[1].values[1,:]));
        end

        for vi=1:length(var)
            components = var[vi].total_components;
            for compi=1:components
                vect[(compi+tmp):totalcomponents:end] = var[vi].values[compi,:];
                tmp = tmp + 1;
            end
        end
    else
        components = var.total_components;
        if vect === nothing
            vect = zeros(components * length(var.values[1,:]));
        end
        for compi=1:components
            vect[compi:components:end] = var.values[compi,:];
        end
    end

    return vect;
end

end #module

ext: jl


text:
import Mimi.compinstance

"""
compute_scc(m::Model=get_model(); year::Int = nothing, gas::Symbol = :CO2, last_year::Int = 3000, equity_weights::Bool = false, eta::Float64 = 1.45, prtp::Float64 = 0.015)

Computes the social cost of CO2 (or other gas if specified) for an emissions pulse in `year`
for the provided MimiFUND model. If no model is provided, the default model from MimiFUND.get_model() is used.
The discount factor is computed from the specified `eta` and pure rate of time preference `prtp`.
Optional regional equity weighting can be used by specifying `equity_weights=true`.
"""
function compute_scc(m::Model=get_model(); year::Union{Int, Nothing} = nothing, gas::Symbol = :CO2, last_year::Int = 3000, equity_weights::Bool = false, eta::Float64 = 1.45, prtp::Float64 = 0.015)

    year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2020)`.") : nothing
    !(last_year in 1950:3000) ? error("Invlaid value of $last_year for last_year. last_year must be within the model's time index 1950:3000.") : nothing
    !(year in 1950:last_year) ? error("Cannot compute the scc for year $year, year must be within the model's time index 1950:$last_year.") : nothing

    mm = get_marginal_model(m; year = year, gas = gas)

    return _compute_scc(mm, year=year, gas=gas, last_year=last_year, equity_weights=equity_weights, eta=eta, prtp=prtp)
end

"""
compute_scc_mm(m::Model=get_model(); year::Union{Int, Nothing} = nothing, gas::Symbol = :CO2, last_year::Int = 3000, equity_weights::Bool = false, eta::Float64 = 1.45, prtp::Float64 = 0.015)

Returns a NamedTuple (scc=scc, mm=mm) of the social cost of carbon and the MarginalModel used to compute it.
Computes the social cost of CO2 (or other gas if specified) for an emissions pulse in `year`
for the provided MimiFUND model. If no model is provided, the default model from MimiFUND.get_model() is used.
The discount factor is computed from the specified `eta` and pure rate of time preference `prtp`.
Optional regional equity weighting can be used by specifying `equity_weights=true`.
"""
function compute_scc_mm(m::Model=get_model(); year::Union{Int, Nothing} = nothing, gas::Symbol = :CO2, last_year::Int = 3000, equity_weights::Bool = false, eta::Float64 = 1.45, prtp::Float64 = 0.015)
    year === nothing ? error("Must specify an emission year. Try `compute_scc_mm(m, year=2020)`.") : nothing
    !(last_year in 1950:3000) ? error("Invlaid value of $last_year for last_year. last_year must be within the model's time index 1950:3000.") : nothing
    !(year in 1950:last_year) ? error("Cannot compute the scc for year $year, year must be within the model's time index 1950:$last_year.") : nothing

    mm = get_marginal_model(m; year = year, gas = gas)
    scc = _compute_scc(mm; year=year, gas=gas, last_year=last_year, equity_weights=equity_weights, eta=eta, prtp=prtp)

    return (scc = scc, mm = mm)
end

# helper function for computing SCC from a MarginalModel, not to be exported
function _compute_scc(mm::MarginalModel; year::Int, gas::Symbol, last_year::Int, equity_weights::Bool, eta::Float64, prtp::Float64)
    ntimesteps = getindexfromyear(last_year)
    run(mm; ntimesteps = ntimesteps)

    marginaldamage = mm[:impactaggregation, :loss]

    ypc = mm.base[:socioeconomic, :ypc]

    # Compute discount factor with or without equityweights
    df = zeros(ntimesteps, 16)
    if !equity_weights
        for r = 1:16
            x = 1.
            for t = getindexfromyear(year):ntimesteps
                df[t, r] = x
                gr = (ypc[t, r] - ypc[t - 1, r]) / ypc[t - 1,r]
                x = x / (1. + prtp + eta * gr)
            end
        end
    else
        globalypc = mm.base[:socioeconomic, :globalypc]
        df = Float64[t >= getindexfromyear(year) ? (globalypc[getindexfromyear(year)] / ypc[t, r]) ^ eta / (1.0 + prtp) ^ (t - getindexfromyear(year)) : 0.0 for t = 1:ntimesteps, r = 1:16]
    end

    # Compute global SCC
    scc = sum(marginaldamage[2:ntimesteps, :] .* df[2:ntimesteps, :]) #* (gas == :C ? 44./12. : 1) # convert to carbon from co2?
    return scc
end

"""
get_marginal_model(m::Model = get_model(); year::Int = nothing, gas::Symbol = :CO2)

Creates a Mimi MarginalModel where the provided m is the base model, and the marginal model has additional emissions of gas `gas` in year `year`.
If no Model m is provided, the default model from MimiFUND.get_model() is used as the base model.
"""
function get_marginal_model(m::Model = get_model(); year::Int = nothing, gas::Symbol = :CO2)
    year == nothing ? error("Must specify emission year. Try `get_marginal_model(m, year=2020)`.") : nothing
    !(year in 1950:3000) ? error("Cannot add marginal emissions in $year, year must be within the model's time index 1950:3000.") : nothing

    mm = create_marginal_model(m, 1e7)  # pulse size is 1MtC for ten years
    add_marginal_emissions!(mm.marginal, year; gas = gas)

    return mm
end

"""
Adds a marginalemission component to m, and sets the additional emissions if a year is specified.
"""
function add_marginal_emissions!(m, year = nothing; gas = :CO2)

    # Add additional emissions to m
    add_comp!(m, Mimi.adder, :marginalemission, before = :climateco2cycle, first = 1951)
    nyears = length(Mimi.time_labels(m))
    addem = zeros(nyears - 1)   # starts one year later, in 1951
    if year != nothing
        addem[getindexfromyear(year)-1:getindexfromyear(year) + 8] .= 1.0
    end
    set_param!(m, :marginalemission, :add, addem)

    # Reconnect the appropriate emissions in m
    if gas == :CO2
        connect_param!(m, :marginalemission, :input, :emissions, :mco2)
        connect_param!(m, :climateco2cycle, :mco2, :marginalemission, :output, repeat([missing], nyears))
    elseif gas == :CH4
        connect_param!(m, :marginalemission, :input, :emissions, :globch4)
        connect_param!(m, :climatech4cycle, :globch4, :marginalemission, :output, repeat([missing], nyears))
    elseif gas == :N2O
        connect_param!(m, :marginalemission, :input, :emissions, :globn2o)
        connect_param!(m, :climaten2ocycle, :globn2o, :marginalemission, :output, repeat([missing], nyears))
    elseif gas == :SF6
        connect_param!(m, :marginalemission, :input, :emissions,:globsf6)
        connect_param!(m, :climatesf6cycle, :globsf6, :marginalemission, :output, repeat([missing], nyears))
    else
        error("Unknown gas: $gas")
    end

end

"""
Helper function to set the marginal emissions in the specified year.
"""
function perturb_marginal_emissions!(m::Model, year; comp_name = :marginalemission)

    ci = compinstance(m, comp_name)
    emissions = Mimi.get_param_value(ci, :add)

    nyears = length(Mimi.dimension(m, :time))
    new_em = zeros(nyears - 1)
    new_em[getindexfromyear(year)-1:getindexfromyear(year) + 8] .= 1.0
    emissions[:] = new_em

end

"""
Creates a MarginalModel of FUND with additional emissions in the specified year for the specified gas.
"""
function create_marginal_FUND_model(; gas = :CO2, year = 2010, parameters = nothing, yearstorun = 1050)

    # Get default FUND model
    FUND = get_model(nsteps = yearstorun, params = parameters)

    # Build marginal model
    mm = create_marginal_model(FUND)
    m1, m2 = mm.base, mm.marginal

    add_marginal_emissions!(m2, year; gas = gas)

    Mimi.build(m1)
    Mimi.build(m2)
    return mm
end

"""
Returns the social cost per one ton of additional emissions of the specified gas in the specified year.
Uses the specified eta and prtp for discounting, with the option to use equity weights.
"""
function get_social_cost(; year = 2010, parameters = nothing, yearstoaggregate = 1000, gas = :CO2, equity_weights = false, eta = 1.0, prtp = 0.001)

    # Get marginal model
    yearstorun = min(1050, getindexfromyear(year) + yearstoaggregate)
    mm = create_marginal_FUND_model(year = year, parameters = parameters, yearstorun = yearstorun, gas = gas)
    run(mm)
    m1, m2 = mm.base, mm.marginal

    damage1 = m1[:impactaggregation, :loss]
    # Take out growth effect effect of run 2 by transforming the damage from run 2 into % of GDP of run 2, and then multiplying that with GDP of run 1
    damage2 = m2[:impactaggregation, :loss] ./ m2[:socioeconomic, :income] .* m1[:socioeconomic, :income]

    # Calculate the marginal damage between run 1 and 2 for each
    # year/region
    marginaldamage = (damage2 .- damage1) / 10000000.0

    ypc = m1[:socioeconomic, :ypc]

    df = zeros(yearstorun + 1, 16)
    if !equity_weights
        for r = 1:16
            x = 1.
            for t = getindexfromyear(year):yearstorun
                df[t, r] = x
                gr = (ypc[t, r] - ypc[t - 1, r]) / ypc[t - 1,r]
                x = x / (1. + prtp + eta * gr)
            end
        end
    else
        globalypc = m1[:socioeconomic, :globalypc]
        df = Float64[t >= getindexfromyear(year) ? (globalypc[getindexfromyear(year)] / ypc[t, r]) ^ eta / (1.0 + prtp) ^ (t - getindexfromyear(year)) : 0.0 for t = 1:yearstorun + 1, r = 1:16]
    end

    scc = sum(marginaldamage[2:end, :] .* df[2:end, :])
    return scc

end

"""
Returns a matrix of marginal damages per one ton of additional emissions of the specified gas in the specified year.
"""
function getmarginaldamages(; year=2010, parameters = nothing, yearstoaggregate = 1000, gas = :CO2)

    # Get marginal model
    yearstorun = min(1050, getindexfromyear(year) + yearstoaggregate)
    mm = create_marginal_FUND_model(year = year, parameters = parameters, yearstorun = yearstorun, gas = gas)
    run(mm)

    # Get damages
    marginaldamages = mm[:impactaggregation, :loss] / 10000000.0
    return marginaldamages
end


text:
## Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
##
## Permission to use, copy, modify, and/or distribute this software for any
## purpose with or without fee is hereby granted, provided that the above
## copyright notice and this permission notice appear in all copies.
##
## THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
## WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
## MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
## ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
## WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
## ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
## OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

struct SoftData
  function SoftData(A, B)
    Base.depwarn("`SoftData(AUX, TI -> AUXTI)` is deprecated, please pass `soft=[(AUX,AUXTI)]` to `iqsim` instead.", :SoftData)
  end
end

const HardData = Dict{NTuple{3,Integer},Real}


text:
make_reader(::Type{PcapBufferReader}, data::Vector{UInt8}) = PcapBufferReader(data)
make_reader(::Type{PcapStreamReader}, data::Vector{UInt8}) = PcapStreamReader(IOBuffer(data))

@testset "$T" for T in (PcapBufferReader, PcapStreamReader)

@testset "empty" begin
    @test_throws EOFError make_reader(T, UInt8[])
end

@testset "invalid" begin
    @test_throws EOFError make_reader(T, [0x0])
    @test_throws ArgumentError make_reader(T, zeros(UInt8, 100))
end

@testset "formats $endianness $timeres" for endianness in (:host, :swapped), timeres in (:us, :ns)
    swapfun = endianness == :swapped ? bswap : identity
    magic = timeres == :ns ? 0xa1b23c4d : 0xa1b2c3d4
    data = zeros(UInt8, 24 + 16 + 8) # file header + record header + payload
    unsafe_store!(Ptr{UInt32}(pointer(data)), swapfun(magic)) # usec
    unsafe_store!(Ptr{UInt16}(pointer(data) + 4), swapfun(UInt16(0x0002))) # ver major
    unsafe_store!(Ptr{UInt16}(pointer(data) + 6), swapfun(UInt16(0x0004))) # ver minor
    unsafe_store!(Ptr{UInt32}(pointer(data) + 16), swapfun(UInt32(65535))) # snaplen
    unsafe_store!(Ptr{UInt32}(pointer(data) + 20), swapfun(LINKTYPE_ETHERNET)) # linktype
    unsafe_store!(Ptr{UInt32}(pointer(data) + 24), swapfun(UInt32(11))) # ts_sec
    unsafe_store!(Ptr{UInt32}(pointer(data) + 28), swapfun(UInt32(13))) # ts_usec
    unsafe_store!(Ptr{UInt32}(pointer(data) + 32), swapfun(UInt32(8))) # incl_len
    unsafe_store!(Ptr{UInt32}(pointer(data) + 36), swapfun(UInt32(8))) # orig_len
    unsafe_store!(Ptr{UInt64}(pointer(data) + 40), 0x0102030405060708) # orig_len
    reader = make_reader(T, data)
    @test reader.header.magic == magic
    @test reader.header.version_major == 2
    @test reader.header.version_minor == 4
    @test reader.header.snaplen == 65535
    @test reader.header.linktype == LINKTYPE_ETHERNET
    record = read(reader)
    @test record.header.ts_sec == 11
    @test record.header.ts_usec == 13
    @test record.header.incl_len == 8
    @test record.header.orig_len == 8
    @test Dates.value(record.timestamp) == 11 * 1_000_000_000 + 13 * (timeres == :ns ? 1 : 1000)
    @test unsafe_load(Ptr{UInt64}(pointer(record.data))) == 0x0102030405060708
    @test eof(reader)
    @test_throws EOFError read(reader)
end

end


text:
# Julia wrapper for header: /usr/include/scip/pub_compr.h
# Automatically generated using Clang.jl wrap_c

function SCIPcomprComp(elem1, elem2)
    ccall((:SCIPcomprComp, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPcomprCompName(elem1, elem2)
    ccall((:SCIPcomprCompName, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPcomprGetData(compr)
    ccall((:SCIPcomprGetData, libscip), Ptr{SCIP_COMPRDATA}, (Ptr{SCIP_COMPR},), compr)
end

function SCIPcomprSetData(compr, comprdata)
    ccall((:SCIPcomprSetData, libscip), Cvoid, (Ptr{SCIP_COMPR}, Ptr{SCIP_COMPRDATA}), compr, comprdata)
end

function SCIPcomprGetName(heur)
    ccall((:SCIPcomprGetName, libscip), Cstring, (Ptr{SCIP_COMPR},), heur)
end

function SCIPcomprGetDesc(compr)
    ccall((:SCIPcomprGetDesc, libscip), Cstring, (Ptr{SCIP_COMPR},), compr)
end

function SCIPcomprGetPriority(compr)
    ccall((:SCIPcomprGetPriority, libscip), Cint, (Ptr{SCIP_COMPR},), compr)
end

function SCIPcomprGetMinNodes(compr)
    ccall((:SCIPcomprGetMinNodes, libscip), Cint, (Ptr{SCIP_COMPR},), compr)
end

function SCIPcomprGetNCalls(compr)
    ccall((:SCIPcomprGetNCalls, libscip), Clonglong, (Ptr{SCIP_COMPR},), compr)
end

function SCIPcomprGetNFound(compr)
    ccall((:SCIPcomprGetNFound, libscip), Clonglong, (Ptr{SCIP_COMPR},), compr)
end

function SCIPcomprIsInitialized(compr)
    ccall((:SCIPcomprIsInitialized, libscip), UInt32, (Ptr{SCIP_COMPR},), compr)
end

function SCIPcomprGetSetupTime(compr)
    ccall((:SCIPcomprGetSetupTime, libscip), Cdouble, (Ptr{SCIP_COMPR},), compr)
end

function SCIPcomprGetTime(compr)
    ccall((:SCIPcomprGetTime, libscip), Cdouble, (Ptr{SCIP_COMPR},), compr)
end


text:
"""
    AbstractFitFunction{T, ND, NP}

Abstract type for an `ND`-dimensional fit of eltype `T` with `NP` parameters.
"""
abstract type AbstractFitFunction{T, ND, NP} end

"""
    struct FitFunction{T, ND, NP} <: AbstractFitFunction{T, ND, NP}

T: Precision type

ND: Dimensionality

NP: NP parameters

Fields:
- `model::Function`: Function of the fit model.
- `fitranges::NTuple{N, Vector{T}}`: Range on which the fit is performed.
- `parameter_names::Vector{Symbol}`: Parameters names.
- `fitted_parameters::Vector{T}`: Fitted parameters.
- `initial_parameters::Vector{T}`: Initial parameters.

"""
struct FitFunction{T, ND, NP} <: AbstractFitFunction{T, ND, NP}
    model::Function
    fitranges::NTuple{ND, Vector{T}}
    parameter_names::Vector{Symbol}
    fitted_parameters::Vector{T}
    initial_parameters::Vector{T}

    function FitFunction{T}(model::Function, ndims::Int, nparams::Int) where {T <: AbstractFloat}
        fitranges::NTuple{ndims, Vector{T}} = NTuple{ndims, Vector{T}}( [-Inf, Inf] for idim in 1:ndims)
        parameter_names::Vector{Symbol} = [Symbol("par$(ipar)") for ipar in 1:nparams]
        fitted_parameters::Vector{T} = [ T(NaN) for ipar in 1:nparams]
        initial_parameters::Vector{T} = [ T(NaN) for ipar in 1:nparams]
        return new{T, ndims, nparams}(model, fitranges, parameter_names, fitted_parameters, initial_parameters)
    end
end

get_ndims(ff::FitFunction{T, ND}) where {T <: AbstractFloat, ND} = ND
get_nparams(ff::FitFunction{T, ND, NP}) where {T <: AbstractFloat, ND, NP} = NP

function set_fitranges!(ff::FitFunction{T, N}, fitranges::NTuple{N, NTuple{2, <:Real}})::Nothing where {T <: AbstractFloat, N}
     for idim in 1:N
        ff.fitranges[idim][:] = T.([ fitranges[idim][1], fitranges[idim][2] ])
    end
    nothing
end
function set_fitranges!(ff::FitFunction{T, N}, fitranges::NTuple{N, AbstractVector{<:Real}})::Nothing where {T <: AbstractFloat, N}
    for idim in 1:N
        @assert length(fitranges[idim]) == 2 "All vectors in `fitranges` must have length 2. But in dimension $(idim) it is $(length(fitranges[idim])) ($(fitranges[idim]))."
        ff.fitranges[idim][:] = T.([ first(fitranges[idim]), last(fitranges[idim]) ])
    end
    nothing
end
function set_parameter_names!(ff::FitFunction{T}, parameter_names::Vector{String})::Nothing where {T <: AbstractFloat}
    ff.parameter_names[:] = [Symbol(name) for name in parameter_names]
    nothing
end
function set_initial_parameters!(ff::FitFunction{T}, initial_parameters::NamedTuple)::Nothing where {T <: AbstractFloat}
    ff.parameter_names[:] = collect(keys(initial_parameters))
    ff.initial_parameters[:] = collect(T, initial_parameters)
    nothing
end
function set_initial_parameters!(ff::FitFunction{T}, initial_parameters::AbstractVector{<:Real})::Nothing where {T <: AbstractFloat}
    ff.initial_parameters[:] = collect(T, initial_parameters)
    nothing
end
function _set_fitted_parameters!(ff::FitFunction{T}, fitted_parameters::AbstractVector{<:Real})::Nothing where {T <: AbstractFloat}
    ff.fitted_parameters[:] = T.(fitted_parameters)
    nothing
end
function _set_fitted_parameters!(ff::FitFunction{T}, fitted_parameters::NamedTuple)::Nothing where {T <: AbstractFloat}
    ff.fitted_parameters[:] = collect(T, fitted_parameters)
    nothing
end

function get_fitted_parameters(ff::FitFunction{T, ND, NP}) where {T <: AbstractFloat, ND, NP}
    return NamedTuple{Tuple(ff.parameter_names), NTuple{NP, T}}( ff.fitted_parameters )
end
function get_initial_parameters(ff::FitFunction{T, ND, NP}) where {T <: AbstractFloat, ND, NP}
    return NamedTuple{Tuple(ff.parameter_names), NTuple{NP, T}}( ff.initial_parameters )
end

function println(io::IO, f::FitFunction)
    println(io, "Model: ", f.model)
    println(io, "  fit ranges: ", f.fitranges)
    println(io, "  Parameters: value (initial value)")
    for i in eachindex(f.parameter_names)
        println("    ", f.parameter_names[i], ": ", f.fitted_parameters[i], " (", f.initial_parameters[i],")" )
    end
end
function print(io::IO, f::FitFunction)
    println(io, f)
end
function show(io::IO, f::FitFunction)
    println(io, f)
end
function show(io::IO, ::MIME"text/plain", f::FitFunction)
    show(io, f)
end

@recipe function f(ff::FitFunction; npoints = 501, use_initial_parameters = false)
    x = collect(range(ff.fitranges[1][1], stop=ff.fitranges[1][2], length=npoints))
    par = use_initial_parameters ? ff.initial_parameters : ff.fitted_parameters
    y = ff.model(x, collect(par))
    linecolor --> (use_initial_parameters ? :green : :red)
    label --> (use_initial_parameters ? "Fit model with initial parameters" : "Fit model with fitted parameters")
    x,y
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct TagCount <: SwaggerModel
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    value::Any # spec type: Union{ Nothing, Int32 } # spec name: value

    function TagCount(;type=nothing, value=nothing)
        o = new()
        validate_property(TagCount, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(TagCount, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type TagCount

const _property_map_TagCount = Dict{Symbol,Symbol}(Symbol("type")=>Symbol("type"), Symbol("value")=>Symbol("value"))
const _property_types_TagCount = Dict{Symbol,String}(Symbol("type")=>"String", Symbol("value")=>"Int32")
Base.propertynames(::Type{ TagCount }) = collect(keys(_property_map_TagCount))
Swagger.property_type(::Type{ TagCount }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_TagCount[name]))}
Swagger.field_name(::Type{ TagCount }, property_name::Symbol) =  _property_map_TagCount[property_name]

function check_required(o::TagCount)
    true
end

function validate_property(::Type{ TagCount }, name::Symbol, val)
end


text:
# using Gadfly
using DataFrames
# using L1AdaptiveTrendFilter

include("../src/CDtypes.jl")
include("../src/initializations.jl")
include("../src/CDmainAlg.jl")
include("../src/GetResult.jl")
include("../src/GramMatrix.jl")
include("../src/InnerProducts.jl")
include("../src/Normalization.jl")
include("../src/utils.jl")

y = readtable(
  "/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/OpticFiber_CaseStudy/ExperimentalData.csv",
  header=false
  )
y = convert(Array, y[:, 2])
y_original = copy(y)

@time BCD, y_path, β_best, y_best, λ_best, γ_best = l1_adaptive_trend_filter(
  y, [1, 3], numλ=100, numγ=3, verbose=true, upper_bounds=[0, 0,0,0,0]
  )

# components
step = β_best[1]
# spike = β_best[2]
slope = β_best[3]

# # plot (Gadfly)
# draw(
#   SVG("/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/OpticFiber_CaseStudy/step.svg", 14inch, 8inch),
#   plot(
#     x=1:length(step), y=step, Geom.point, Geom.line,
#     Guide.xlabel("Hours"), Guide.ylabel("Step size"),
#     Coord.Cartesian(xmin=0,xmax=length(y_original))
#     ))
# # draw(
# #   SVG("/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/OpticFiber_CaseStudy/spike.svg", 14inch, 8inch),
# #   plot(x=1:length(spike),y=spike, Geom.point, Geom.line, Guide.xlabel("Hours")))
# draw(
#   SVG("/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/OpticFiber_CaseStudy/slope.svg", 14inch, 8inch),
#   plot(x=1:length(slope),y=slope, Geom.point, Geom.line)
#   )
# draw(
#   SVG("/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/OpticFiber_CaseStudy/fit.svg", 7inch, 4.5inch),
#   plot(
#     layer(x=1:length(y_original), y=y_original, Geom.line, Theme(default_color=color("red"))),
#     layer(x=1:length(y_best), y=y_best, Geom.line, Theme(default_color=color("blue"))),
#     Guide.xlabel("Hours"), Guide.ylabel("MW"), Guide.title("Wind Power Fit"),
#     Coord.Cartesian(xmin=0,xmax=8000)
#     ))

# write results
writecsv("/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/OpticFiber_CaseStudy/y_best.csv", y_best)
writecsv("/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/OpticFiber_CaseStudy/step.csv", step)

max_stars_repo_name: joaquimg/ComponentFilterCD

max_issues_repo_name: joaquimg/ComponentFilterCD

max_forks_repo_name: joaquimg/ComponentFilterCD


text:
using BinaryBuilder, Pkg

name = "SCALAPACK32"
version = v"2.1.0"

sources = [
  ArchiveSource("http://www.netlib.org/scalapack/scalapack-$(version).tgz",
                "61d9216cf81d246944720cfce96255878a3f85dec13b9351f1fa0fd6768220a6"),
  DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
mkdir -p ${libdir}
cd $WORKSPACE/srcdir/scalapack-*

# the patch prevents running foreign executables, which fails on most platforms
# we instead set CDEFS manually below
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
  atomic_patch -p1 ${f}
done

CMAKE_FLAGS=(-DCMAKE_INSTALL_PREFIX=${prefix} \
             -DCMAKE_FIND_ROOT_PATH=${prefix} \
             -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
             -DCMAKE_BUILD_TYPE=Release \
             -DBLAS_LIBRARIES="-lopenblas" \
             -DLAPACK_LIBRARIES="-lopenblas" \
             -DBUILD_SHARED_LIBS=ON \
             -DMPIEXEC="${bindir}/mpirun")

if [[ "${target}" == i686-*  ]] || [[ "${target}" == x86_64-*  ]]; then
  CMAKE_FLAGS+=(-DCMAKE_EXE_LINKER_FLAGS="-lgfortran -lquadmath")
else
  CMAKE_FLAGS+=(-DCMAKE_EXE_LINKER_FLAGS="-lgfortran")
fi

if [[ "${target}" == aarch64-apple-darwin* ]]; then
  CMAKE_FLAGS+=(-DCMAKE_Fortran_FLAGS="-fallow-argument-mismatch")
fi

export CDEFS="Add_"

mkdir build
cd build
cmake .. "${CMAKE_FLAGS[@]}"

make -j${nproc} all
make install

mv -v ${libdir}/libscalapack.${dlext} ${libdir}/libscalapack32.${dlext}

PATCHELF_FLAGS=()

# ppc64le and aarch64 have 64KB page sizes, don't muck up the ELF section load alignment
if [[ ${target} == aarch64-* || ${target} == powerpc64le-* ]]; then
  PATCHELF_FLAGS+=(--page-size 65536)
fi

if [[ ${target} == *linux* ]] || [[ ${target} == *freebsd* ]]; then
  patchelf ${PATCHELF_FLAGS[@]} --set-soname libscalapack32.${dlext} ${libdir}/libscalapack32.${dlext}
elif [[ ${target} == *apple* ]]; then
  install_name_tool -id libscalapack32.${dlext} ${libdir}/libscalapack32.${dlext}
fi
"""

# OpenMPI and MPICH are not precompiled for Windows
# Can't get the code to build for PowerPC with libgfortran3
platforms = expand_gfortran_versions(supported_platforms(; exclude=p -> Sys.iswindows(p) || arch(p) == "powerpc64le"))

# The products that we will ensure are always built
products = [
    LibraryProduct("libscalapack32", :libscalapack32),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="OpenBLAS32_jll", uuid="656ef2d0-ae68-5445-9ca0-591084a874a2")),
    Dependency(PackageSpec(name="MPICH_jll", uuid="7cb0a576-ebde-5e09-9194-50597f1243b4")),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies, julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    qstat(server)

If sbatch is running on server it will return the output of the `qstat` command.
"""
qstat(server) = server=="localhost" ? run(`qstat`) : sshcmd(server, "qstat")
qstat()       = qstat(getdefault_server())

"""
    watch_qstat(server)
If sbatch is running on server it will return the output of the `watch qstat` command.
"""
watch_qstat(server) = server=="localhost" ? run(`watch qstat`) : sshcmd(server, "watch qstat")
watch_qstat()       = watch_qstat(getdefault_server())

#--------------- Slurm Interactions -------------------#
"""
    slurm_history_jobdir(stardate=yesterday())

Returns the unique job directories of the jobs that ran since the `startdate`.
Startdate should be printed in following format: yyyy-mm-dd.
"""
function slurm_history_jobdir(startdate=yesterday()) #format of startdate = yyyy-mm-dd
    history = slurm_process_command(`sacct --starttime $startdate --format=Workdir%100`)
    output = String[]
    for h in history
        if h ∉ output && ispath(h)
            push!(output, h)
        end
    end
    return output
end

"""
    slurm_jobid(job::DFJob, startdate=yesterday())

Looks through the jobs since the `startdate` and returns the job ID if found.
Returns -1 if the jobID was not found in the list of jobs since `startdate`.
"""
function slurm_jobid(job::DFJob, startdate=yesterday())
    runslocal_assert(job)
    if haskey(job.metadata, :slurmid)
        return job.metadata[:slurmid]
    end
    id_dir = filter(x->length(x) ==2, split.(slurm_process_command(`sacct --starttime $startdate --format=JobID,Workdir%100`)))
    id_ = -1
    for (id, dir) in id_dir
        if dir == job.local_dir[1:end-1]
            id_ = parse(Int, id)
            break
        end
    end
    if id_ == -1
        @info "Job in directory $(job.local_dir) was not found in the slurm jobs since
        $startdate"
    else
        job.metadata[:slurmid] = id_
    end
    return id_
end

"""
    slurm_isrunning(job::DFJob)

Returns whether the job is running.
"""
function slurm_isrunning(job::DFJob)
    runslocal_assert(job)
    id = slurm_jobid(job)
    if id != -1
        if slurm_process_command(`sacct -j $id --format=state`)[1] == "RUNNING"
            return true
        else
            return false
        end
    else
        return false
    end
end

"""
    slurm_isqueued(job::DFJob)

Returns whether the job is queued.
"""
function slurm_isqueued(job::DFJob)
    runslocal_assert(job)
    id = slurm_jobid(job)
    if id != -1
        if slurm_process_command(`sacct -j $id --format=state`)[1] == "PENDING"
            return true
        else
            return false
        end
    else
        return false
    end
end

"""
    slurm_mostrecent(index=1, jobfile="job.tt", startdate=lastmonth(), args...; kwargs...)

Returns whether the `index`th most recent job with job script `jobfile`.
Extra args and kwargs will be passed to the `DFJob` constructor.
"""
function slurm_mostrecent(index=1, jobfile="job.tt", startdate=lastmonth(), args...; kwargs...)
    dirs = slurm_history_jobdir(startdate)
    @assert length(dirs) >= index "There are less recent jobs since startdate $startdate than the index $index."
    jobpath = joinpath(dirs[index], jobfile)
    @assert ispath(jobpath) "The directory $(dirs[index]) does not have a job file with filename $jobfile."
    return DFJob(dirs[index], args...; job_fuzzy=jobfile, kwargs...)
end


text:
struct binned1dDensity{T}
    grid::Vector{Float64}
    cumarr::Vector{T}
    density::Function
end

function binned1dDensity(grid, weights; density::Function = x->(println("warning: no callback provided!");1))
    _weights = weights ./ sum(weights)
    cumarr = [0;cumsum(vcat(_weights...), dims=1)]
    return binned1dDensity(grid, cumarr, density)
end

function getbinned1dDensity(g, lims, Nbins)
    grid = collect(range(lims[1], lims[2], length=Nbins))
    #
    weights = [g(v) for v in (grid[2:end] .+ grid[1:end-1]) ./ 2]
    return binned1dDensity(grid, weights; density=g)
end

import Base:rand
function rand(bD::binned1dDensity; removenegative::Bool=false)
    # get
    binind = findfirst(bD.cumarr .> rand())-1
    σl, σr = bD.grid[binind], bD.grid[binind+1]
    σ = σl+rand()*(σr-σl)
    removenegative && bD.density(σ) == 0.0 && return rand(bD,removenegative)
    return σ
end

function generate(d::T where T <: AbstractFunctionWithParameters, N::Int;
    p=freepars(d), Nbins=100, removenegative::Bool=false)
    grid = collect(range(lims(d)..., length=Nbins))
    centers = (grid[2:end] .+ grid[1:end-1]) ./ 2
    #
    weights = d(centers; p)
    bD = binned1dDensity(grid, weights; density=x->d(x;p))
    return [rand(bD; removenegative) for _ in 1:N]
end
rand(d::T where T <: AbstractFunctionWithParameters, N::Int=1) = generate(d, N)


text:
# Examples
"""

R=2
C=3

println("2D pseudo-index fill ...")
A=init_nested_blocks(R,C)
for c=1:C,r=1:R
    put_nested_block!(randn(1,1),A,r,c)
    # or just A[r][c]=r*c
end
display(A);println()
display(A[1]);println()
display(A[1][1]);println()

println("linear fill column-wise ...")
A=init_nested_blocks(R,C)
for i=1:R*C put_nested_block!(string(i),A,i) end
display(A);println()
display(A[1]);println()
display(A[1][1]);println()

println("linear fill row-wise")
A=init_nested_blocks(R,C)
for i=1:R*C put_nested_block!(string(i),A,i,colwise=false) end
display(A);println()
display(A[1]);println()
display(A[1][1]);println()

"""

export init_nested_blocks
"""
    julia> A = init_nested_blocks(R,C)

Initialize empty Vector{Vector} - pseudo-2D nested block array with

- R number of rows
- C number of columns

Row blocks can be accessed using A[r]

Single blocks can be accessed with A[r][c]

"""
function init_nested_blocks(R::Integer,C::Integer)
    A=Vector{Vector}(undef,R)
    for r=1:R
        A[r]=Vector{Any}(undef,C)
    end
    return A
end

export put_nested_block!
"""
    julia> put_nested_block!(O,A,r,c)

Assign any object to element of Vector{Vector} created with `init_nested_blocks`, where

- O is any object
- A is Vector{Vector} from `init_nested_blocks`
- r is row index
- c is column index

This is equivalent to A[r][c]=O

"""
function put_nested_block!(O,A::Vector{Vector},r::Integer,c::Integer)
    A[r][c]=O
    return nothing
end

export put_nested_block_row_wise!
"""
    julia> put_nested_block_row_wise!(O,A,l)

Assign any object to element of Vector{Vector} created with `init_nested_blocks`, uing linear index in row-wise convention, where

- O is any object
- A is Vector{Vector} from `init_nested_blocks`
- l is linear index in row-wise convention

"""
function put_nested_block_row_wise!(O,A,l::Integer)
    R=size(A,1)
    C=size(A[1],1)
    @assert l > 0
    @assert l <= R*C
    r=ceil(Int,l/C)
    c=(l-1)%C+1
    A[r][c]=O
    return nothing
end

export put_nested_block_col_wise!
"""
    julia> put_nested_block_col_wise!(O,A,l)

Assign any object to element of Vector{Vector} created with `init_nested_blocks`, uing linear index in column-wise convention, where

- O is any object
- A is Vector{Vector} from `init_nested_blocks`
- l is linear index in column-wise convention

"""
function put_nested_block_col_wise!(O,A::Vector{Vector},l::Integer)
    R=size(A,1)
    C=size(A[1],1)
    @assert l > 0
    @assert l <= R*C
    c=ceil(Int,l/R)
    r=(l-1)%R+1
    A[r][c]=O
    return nothing
end

export put_nested_block!
"""
    julia> put_nested_block!(O,A,l;colwise=true)

Assign any object to element of Vector{Vector} created with `init_nested_blocks`, uing linear index in selected convention, where

- O is any object
- A is Vector{Vector} from `init_nested_blocks`
- l is linear index in selected convention
- colwise keyword
    - true uses column-wise convention (default)
    - false uses row-wise convention

"""
function put_nested_block!(O,A::Vector{Vector},l::Integer;colwise::Bool=true)
    if colwise
        put_nested_block_col_wise!(O,A,l)
    else
        put_nested_block_row_wise!(O,A,l)
    end
end


text:
struct ApproxFunProblem{P,S} <: DEProblem
    prob::P
    space::S
end

function ApproxFunProblem(prob;n = ncoefficients(prob.u0))
  sp = space(prob.u0)
  _prob = ODEProblem((u,p,t)->pad!(prob.f(Fun(sp,u),p,t).coefficients,n),pad!(prob.u0.coefficients,n),prob.tspan,prob.p)
  return ApproxFunProblem(_prob,sp)
end

struct ApproxFunSolution{SOL,S}
    sol::SOL
    space::S
end

(sol::ApproxFunSolution)(t) = Fun(sol.space,sol.sol(t))
(sol::ApproxFunSolution)(t,x) = sol(t)(x)

@inline Base.getindex(sol::ApproxFunSolution, I::Int...) = Fun(sol.space,sol.sol[I[end]])[Base.front(I)...]
@inline Base.getindex(sol::ApproxFunSolution, I::Int) = Fun(sol.space,sol.sol[I])
@inline Base.getindex(sol::ApproxFunSolution, I::Colon) = Fun.(sol.space,sol.sol[I])
@inline Base.getindex(sol::ApproxFunSolution, I::AbstractArray{Int}) = Fun.(sol.space,sol.sol[I])

@inline Base.length(sol::ApproxFunSolution) = length(sol.sol)
@inline function Base.endof(sol::ApproxFunSolution)
    endof(sol.sol)
end
@inline Base.size(sol::ApproxFunSolution) = size(sol.sol)
@inline Base.eachindex(sol::ApproxFunSolution) = Base.OneTo(length(sol.u))

function DiffEqBase.init(prob::ApproxFunProblem,alg,args...;kwargs...)
    init(prob.prob,alg,args...;kwargs...)
end
function DiffEqBase.solve(prob::ApproxFunProblem,alg,args...;kwargs...)
    sol = solve(prob.prob,alg,args...;kwargs...)
    ApproxFunSolution(sol,prob.space)
end

function Base.show(io::IO, A::ApproxFunProblem)
  print(io,"Problem: ")
  show(io,A.prob)
  println(io)
  print(io,"Space: ")
  show(io, A.space)
end


text:
println("Hello World")
hello world()
5y4gbnbvcvbnhgfdu vfnpf

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
function MA.mutable_operate!(op::MA.AddSubMul, aff::PAE, x::Union{JuMP.VariableRef, GAEv, ParameterRef, GAEp}, c::Number)
    return MA.mutable_operate!(op, aff, c, x)
end
function MA.mutable_operate!(op::MA.AddSubMul, aff::PAE, c::Number, x::Union{JuMP.VariableRef, GAEv})
    if !iszero(c)
        MA.mutable_operate!(op, aff.v, c, x)
    end
    return aff
end
function MA.mutable_operate!(op::MA.AddSubMul, aff::PAE, c::Number, x::Union{ParameterRef, GAEp})
    if !iszero(c)
        MA.mutable_operate!(op, aff.p, c, x)
    end
    return aff
end
function MA.mutable_operate!(op::MA.AddSubMul, aff::PAE, c::Number, x::Number)
    if !iszero(c) && !iszero(x)
        aff.v = MA.mutable_operate!(op, aff.v, c, x)
    end
    return aff
end

function JuMP.add_to_expression!(aff::PAE, other::Number)
    JuMP.add_to_expression!(aff.v, other)
end
function JuMP.add_to_expression!(aff::PAE, new_var::JuMP.VariableRef, new_coef)
    JuMP.add_to_expression!(aff.v, new_coef, new_var)
end
function JuMP.add_to_expression!(aff::PAE, new_coef, new_var::JuMP.VariableRef)
    JuMP.add_to_expression!(aff.v, new_coef, new_var)
end
function JuMP.add_to_expression!(aff::PAE, new_var::JuMP.VariableRef)
    JuMP.add_to_expression!(aff.v, new_var)
end
function JuMP.add_to_expression!(aff::PAE, new_param::ParameterRef, new_coef)
    JuMP.add_to_expression!(aff.p, new_coef, new_param)
end
function JuMP.add_to_expression!(aff::PAE, new_param::ParameterRef)
    JuMP.add_to_expression!(aff.p,new_param)
end
function JuMP.add_to_expression!(aff::PAE, new_coef, new_param::ParameterRef)
    JuMP.add_to_expression!(aff.p, new_coef, new_param)
end
function JuMP.add_to_expression!(lhs_aff::PAE, rhs_aff::PAE)
    JuMP.add_to_expression!(lhs_aff.p, rhs_aff.p)
    JuMP.add_to_expression!(lhs_aff.v, rhs_aff.v)
end


text:
"""
A simple null model consisting of a fixed size population of
organisms where in each generation, some random mutations are
introduced.
"""
struct RandomModel
    current_population::Vector{Vector{Bool}}

    genotype_length::UInt8
    mutation_probability::FLoat64
    population_size::UInt8

    function RandomModel(genotype_length, mutation_probability, population_size)
        new([falses(genotype_length) for _ in 1:population_size],
            genotype_length, mutation_probability, population_size)
    end
end

function execute!(model::RandomModel)
    for organism in model.current_population
        if rand() < model.mutation_probability
            index_to_change = rand(1:model.population_size)
            organism[index_to_change] = !organism[index_to_change]
        end
    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
# prints the classification result produce by running mil_data.jl
# as an input, it takes the folder with saved outputs and possibly a number
# of an output file that only needs to be read
# run from the command line, e.g.
# julia path/to/file/read_output.jl path/with/outputs - this reads all the files in output path
# or
# julia path/to/file/read_output.jl path/with/outputs 1 3 7- this reads first, third and seventh file in the folder
using JLD
include("mil_util.jl")

output_path = ARGS[1]
files = readdir(output_path)
nfiles = size(files)[1]

if size(ARGS)[1] > 1
    inds = map(x->parse(Float64,x),ARGS[2:end])
else
    inds = 1:nfiles
end

for n in inds
    n = convert(Int, n)
    file = files[n]
    println("File n = $n")
    try
        data = load("$output_path/$file")
        table_summary(data)
    catch e
        println("Something bad happened, probably wrong file type.")
        println(e)
    end
    n+=1
end


text:
module Datasets

using DataDeps
using HTTP
using InternedStrings
using WordTokenizers

const Container{T} = Union{NTuple{N, T}, Vector{T}} where N

export Dataset, Train, Dev, Test
export dataset, datafile, get_batch, get_vocab, get_labels, batched

include("download_utils.jl")
include("./dataset.jl")

include("translate/wmt.jl")
using .WMT

include("translate/iwslt.jl")
using .IWSLT

include("qa/clozetest.jl")
using .ClozeTest

include("glue/glue.jl")
using .GLUE

end


text:
if VERSION > v"0.5-" && isdefined(Base,:displaysize)
    tty_cols = Base.displaysize()[2]
else
    tty_cols = Base.tty_size()[2]
end

function set_tty_cols(val::Integer)
  global tty_cols = val
end

show_backtraces = true

function set_show_backtraces(val::Bool)
  global show_backtraces = val
end

function stringify_error(err, bt)
  const buff = IOBuffer()
  if show_backtraces
    showerror(buff, err, bt)
  else
    showerror(buff, err)
  end
  return UTF8String(buff.data)
end

function centre(message::AbstractString, padwith::Char='=')
  const padding = max(0, tty_cols - (length(message)+2))
  const padstr = "$padwith"
  return "$(padstr^round(Int,ceil(padding/2))) $message $(padstr^round(Int,floor(padding/2)))"
end

function underline(message::AbstractString, underlinewith::Char='=')
  const underlinestr = "$underlinewith"
  return "$message\n$(underlinestr^length(message))"
end

function breakwater(padwith::Char='=')
  return "$padwith"^tty_cols
end


text:
function generate_execute(domain::Domain, state::State,
                          domain_type::Symbol, state_type::Symbol)
    execute_def = quote
        function execute(domain::$domain_type, state::$state_type, term::Term)
            execute(domain, state, get_action(domain, term.name), term.args)
        end
        function execute(domain::$domain_type, state::$state_type, term::Term;
                         check::Bool=false)
            execute(domain, state, get_action(domain, term.name), term.args;
                    check=check)
        end
        function execute!(domain::$domain_type, state::$state_type, term::Term)
            execute!(domain, state, get_action(domain, term.name), term.args)
        end
        function execute!(domain::$domain_type, state::$state_type, term::Term;
                          check::Bool=false)
            execute!(domain, state, get_action(domain, term.name), term.args;
                     check=check)
        end
    end
    return execute_def
end

function generate_execute(domain::Domain, state::State,
                          domain_type::Symbol, state_type::Symbol,
                          action_name::Symbol, action_type::Symbol)
    action = get_actions(domain)[action_name]
    varmap = Dict{Var,Any}(a => :(args[$i].name) for (i, a) in
                           enumerate(get_argvars(action)))
    # TODO: Fix mutating execute
    precond = generate_check_expr(domain, state, get_precond(action), varmap,
                                  :prev_state)
    effect = generate_effect_expr(domain, state, get_effect(action), varmap)
    effect! = generate_effect_expr(domain, state, get_effect(action), varmap,
                                   :prev_state)
    execute_def = quote
        function execute(domain::$domain_type, prev_state::$state_type,
                         action::$action_type, args; check::Bool=false)
            if check && !(@inbounds $precond)
                error("Precondition not satisfied")
            end
            state = copy(prev_state)
            @inbounds $effect
            return state
        end
        function execute!(domain::$domain_type, prev_state::$state_type,
                          action::$action_type, args; check::Bool=false)
            if check && !(@inbounds $precond)
                error("Precondition not satisfied")
            end
            @inbounds $effect!
            return prev_state
        end
    end
    return execute_def
end


text:
using Documenter, GXBeam

makedocs(;
    modules = [GXBeam],
    pages = [
        "Home" => "index.md",
        "Getting Started" => "guide.md",
        "Examples" => "examples.md",
        "Using GXBeam with DifferentialEquations.jl" => "diffeq.md",
        "Library" => "library.md"
    ],
    sitename = "GXBeam.jl",
    authors = "Taylor McDonnell <taylormcd@byu.edu>",
)

deploydocs(
    repo = "github.com/byuflowlab/GXBeam.jl.git",
)


text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля, на котором могут находиться также внутренние
#прямоугольные перегородки (все перегородки изолированы друг от друга, прямоугольники могут вырождаться в отрезки)
#РЕЗУЛЬТАТ: Робот - в исходном положении и в углах поля стоят маркеры
include("library.jl")
function putmarkerInCorner!(r::Robot)
    for side in instances(HorizonSide)
        putmarker!(r)
        while !isborder(r, side)
            move!(r, side)
        end
    end
end
function ugly_with_obstacles(r::Robot)
    sides = moveAndReturnDirections!(r)
    putmarkerInCorner!(r)
    moveToStartplace!(r)
    moveToBeginplace!(r, sides)

end

max_stars_repo_name: listova-nastya/mirea-progs

max_issues_repo_name: listova-nastya/mirea-progs

max_forks_repo_name: listova-nastya/mirea-progs


text:
import Base.<, Base.>, Base.-

export PieceColor, PieceType, Piece

export colorfromchar,
    coloropp,
    isok,
    isslider,
    pcolor,
    piecefromchar,
    piecetypefromchar,
    ptype,
    tochar,
    tounicode
export BISHOP,
    BLACK,
    COLOR_NONE,
    EMPTY,
    KING,
    KNIGHT,
    PAWN,
    PIECE_BB,
    PIECE_BK,
    PIECE_BN,
    PIECE_BP,
    PIECE_BQ,
    PIECE_BR,
    PIECE_TYPE_NONE,
    PIECE_WB,
    PIECE_WK,
    PIECE_WN,
    PIECE_WP,
    PIECE_WQ,
    PIECE_WR,
    QUEEN,
    ROOK,
    WHITE

"""
    PieceColor

Type representing the color of a chess piece.

The possible values are `WHITE`, `BLACK` and `COLOR_NONE`. The reason for the
existence of the value `COLOR_NONE` is that we represent a chess board as an
array of pieces, and we need a special `Piece` value `EMPTY` to indicate an
empty square on the board. The color of the `EMPTY` piece is `COLOR_NONE`.
"""
struct PieceColor
    val::Int
end

const WHITE = PieceColor(1)
const BLACK = PieceColor(2)
const COLOR_NONE = PieceColor(3)

function Base.show(io::IO, c::PieceColor)
    if c == WHITE
        print(io, "WHITE")
    elseif c == BLACK
        print(io, "BLACK")
    elseif c == COLOR_NONE
        print(io, "COLOR_NONE")
    else
        print(io, "PieceColor(???)")
    end
end

"""
    isok(c::PieceColor)

Tests whether a `PieceColor` value is a valid.

Returns `true` if and only if `c` is either `WHITE` or `BLACK`.

# Examples

```julia-repl
julia> isok(WHITE)
true

julia> isok(BLACK)
true

julia> isok(COLOR_NONE)
false

julia> isok(PieceColor(42))
false
```
"""
function isok(c::PieceColor)::Bool
    c == WHITE || c == BLACK
end

"""
    coloropp(c::PieceColor)

Returns the opposite of a color.

# Examples

```julia-repl
julia> coloropp(WHITE) == BLACK
true

julia> coloropp(BLACK) == WHITE
true
```
"""
function coloropp(c::PieceColor)::PieceColor
    PieceColor(c.val ⊻ 3)
end

"""
    -(c::PieceColor)

Returns the opposite of a color. Synonymous to `coloropp`.

# Examples
```julia-repl
julia> -WHITE
BLACK

julia> -BLACK
WHITE
```
"""
function -(c::PieceColor)
    coloropp(c)
end

"""
    colorfromchar(c::Char)

Tries to convert a character to a `PieceColor`.

The return value is a `Union{PieceColor, Nothing}`. If the input character is
one of the four characters `'w'`, `'b'`, `'W'`, `'B'`, the function returns the
obvious corresponding color (`WHITE` or `BLACK`). For all other input
characters, the function returns `nothing`.

# Examples

```julia-repl
julia> colorfromchar('w') == WHITE
true

julia> colorfromchar('B') == BLACK
true

julia> colorfromchar('x') == nothing
true
```
"""
function colorfromchar(c::Char)::Union{PieceColor,Nothing}
    i = findfirst(isequal(lowercase(c)), "wb")
    if isnothing(i)
        nothing
    else
        PieceColor(i)
    end
end

"""
    tochar(c::PieceColor)

Converts a color to a character.

# Examples

```julia-repl
julia> tochar(WHITE)
'w': ASCII/Unicode U+0077 (category Ll: Letter, lowercase)

julia> tochar(BLACK)
'b': ASCII/Unicode U+0062 (category Ll: Letter, lowercase)

julia> tochar(COLOR_NONE)
'?': ASCII/Unicode U+003f (category Po: Punctuation, other)
```
"""
function tochar(c::PieceColor)::Char
    if c == WHITE
        'w'
    elseif c == BLACK
        'b'
    else
        '?'
    end
end

"""
    PieceType

Type representing the type of a chess piece.

This is essentially a piece without color. The possible values are `PAWN`,
`KNIGHT`, `BISHOP`, `ROOK`, `QUEEN`, `KING` and `PIECE_TYPE_NONE`. The reason
for the existence of the value `PIECE_TYPE_NONE` is that we represent a chess
board as an array of pieces, and we need a special `Piece` value `EMPTY` to
indicate an empty square on the board. The type of the `EMPTY` piece is
`PIECE_TYPE_NONE`
"""
struct PieceType
    val::Int
end

(<)(t1::PieceType, t2::PieceType) = t1.val < t2.val

const PAWN = PieceType(1)
const KNIGHT = PieceType(2)
const BISHOP = PieceType(3)
const ROOK = PieceType(4)
const QUEEN = PieceType(5)
const KING = PieceType(6)
const PIECE_TYPE_NONE = PieceType(7)

function Base.show(io::IO, t::PieceType)
    if t == PAWN
        print(io, "PAWN")
    elseif t == KNIGHT
        print(io, "KNIGHT")
    elseif t == BISHOP
        print(io, "BISHOP")
    elseif t == ROOK
        print(io, "ROOK")
    elseif t == QUEEN
        print(io, "QUEEN")
    elseif t == KING
        print(io, "KING")
    elseif t == PIECE_TYPE_NONE
        print(io, "PIECE_TYPE_NONE")
    else
        print(io, "PieceType($(t.val))")
    end
end

"""
    isok(t::PieceType)

Tests whether a `PieceType` value is valid.

Returns `true` for all of `PAWN`, `KNIGHT`, `BISHOP`, `ROOK`, `QUEEN` and
`KING`, and `false` for all other inputs.

# Examples

```julia-repl
julia> isok(QUEEN)
true

julia> isok(KNIGHT)
true

julia> isok(PIECE_TYPE_NONE)
false

julia> isok(PieceType(-1))
false
```
"""
function isok(t::PieceType)::Bool
    t >= PAWN && t <= KING
end

"""
    piecetypefromchar(c::Chars)

Tries to convert a character to a `PieceType`.

The return value is a `Union{PieceType, Nothing}`. If the input character is a
valid upper- or lowercase English piece letter (PNBRQK), the function returns
the corresponding piece type. For all other input characters, the function
returns `nothing`.

# Examples

```julia-repl
julia> piecetypefromchar('n') == KNIGHT
true

julia> piecetypefromchar('B') == BISHOP
true

julia> piecetypefromchar('a') == nothing
true
```
"""
function piecetypefromchar(c::Char)::Union{PieceType,Nothing}
    i = findfirst(isequal(lowercase(c)), "pnbrqk")
    if isnothing(i)
        nothing
    else
        PieceType(i)
    end
end

"""
    tochar(t::PieceType, uppercase = false)

Converts a `PieceType` value to a character.

A valid piece type value is converted to its standard English algebraic
notation piece letter. Any invalid piece type value is converted to a `'?'`
character. The optional parameter `uppercase` controls whether the character
is an upper- or lower-case letter.

# Examples

```julia-repl
julia> tochar(PAWN)
'p': ASCII/Unicode U+0070 (category Ll: Letter, lowercase)

julia> tochar(ROOK, true)
'R': ASCII/Unicode U+0052 (category Lu: Letter, uppercase)

julia> tochar(PIECE_TYPE_NONE)
'?': ASCII/Unicode U+003f (category Po: Punctuation, other)
```
"""
function tochar(t::PieceType, uppercase = false)::Char
    if isok(t)
        (uppercase ? "PNBRQK" : "pnbrqk")[t.val]
    else
        '?'
    end
end

"""
    Piece

Type representing a chess piece.

The possible values are `PIECE_WP`, `PIECE_WN`, `PIECE_WB`, `PIECE_WR`,
`PIECE_WQ`, `PIECE_WK`, `PIECE_BP`, `PIECE_BN`, `PIECE_BB`, `PIECE_BR`,
`PIECE_BQ`, `PIECE_BK` and `EMPTY`. The reason for the existence of the
value `EMPTY` is that we represent a chess board as an array of pieces, and
we need a value to indicate an empty square on the board.
"""
struct Piece
    val::Int
end

"""
    Piece(c::PieceColor, t::PieceType)

Construct a piece with the given color and type.

# Examples
```julia-repl
julia> Piece(BLACK, QUEEN)
PIECE_BQ
```
"""
Piece(c::PieceColor, t::PieceType) = Piece(((c.val - 1) << 3) | t.val)

const PIECE_WP = Piece(WHITE, PAWN)
const PIECE_WN = Piece(WHITE, KNIGHT)
const PIECE_WB = Piece(WHITE, BISHOP)
const PIECE_WR = Piece(WHITE, ROOK)
const PIECE_WQ = Piece(WHITE, QUEEN)
const PIECE_WK = Piece(WHITE, KING)
const PIECE_BP = Piece(BLACK, PAWN)
const PIECE_BN = Piece(BLACK, KNIGHT)
const PIECE_BB = Piece(BLACK, BISHOP)
const PIECE_BR = Piece(BLACK, ROOK)
const PIECE_BQ = Piece(BLACK, QUEEN)
const PIECE_BK = Piece(BLACK, KING)
const EMPTY = Piece(COLOR_NONE, PIECE_TYPE_NONE)

function Base.show(io::IO, p::Piece)
    if isok(p)
        print(io, "PIECE_$(uppercase(tochar(pcolor(p))))$(uppercase(tochar(ptype(p))))")
    elseif p == EMPTY
        print(io, "EMPTY")
    else
        print(io, "PieceType($(p.val))")
    end
end

"""
    pcolor(p::Piece)

Find the color of a `Piece`.

# Examples

```julia-repl
julia> pcolor(PIECE_WB)
WHITE

julia> pcolor(EMPTY)
COLOR_NONE
```
"""
function pcolor(p::Piece)::PieceColor
    PieceColor(p.val >> 3 + 1)
end

"""
    ptype(p::Piece)

Find the type of a `Piece`.

# Examples

```julia-repl
julia> ptype(PIECE_BQ)
QUEEN

julia> ptype(EMPTY)
PIECE_TYPE_NONE
```
"""
function ptype(p::Piece)::PieceType
    PieceType(p.val & 7)
end

"""
    isok(p::Piece)

Tests wheter a `Piece` value is valid.

Returns `true` if and only if the color and the type of the piece are valid
piece color and piece type values, respectively.

# Examples

```julia-repl
julia> isok(PIECE_WB)
true

julia> isok(PIECE_BQ)
true

julia> isok(EMPTY)
false

julia> isok(Piece(-10))
false
```
"""
function isok(p::Piece)::Bool
    isok(ptype(p)) && isok(pcolor(p))
end

"""
    piecefromchar(ch::Char)

Tries to convert a character to a `Piece`.

The return value is a `Union{Piece, Nothing}`. If the input character is a
valid English piece letter, the corresponding piece is returned. If the piece
letter is uppercase, the piece is white. If the piece letter is lowercase,
the piece is black.

If the input value is not a valid English piece letter, the function returns
`nothing`.

# Examples

```julia-repl
julia> piecefromchar('Q')
PIECE_WQ

julia> piecefromchar('n')
PIECE_BN

julia> piecefromchar('-') == nothing
true
```
"""
function piecefromchar(ch::Char)::Union{Piece,Nothing}
    c = isuppercase(ch) ? WHITE : BLACK
    t = piecetypefromchar(ch)
    if isnothing(t)
        nothing
    else
        Piece(c, t)
    end
end

"""
    tochar(p::Piece)

Converts a piece to a character.

# Examples

```julia-repl
julia> tochar(PIECE_WN)
'N': ASCII/Unicode U+004e (category Lu: Letter, uppercase)

julia> tochar(PIECE_BK)
'k': ASCII/Unicode U+006b (category Ll: Letter, lowercase)

julia> tochar(EMPTY)
'?': ASCII/Unicode U+003f (category Po: Punctuation, other)
```
"""
function tochar(p::Piece)::Char
    if isok(p)
        tochar(ptype(p), pcolor(p) == WHITE)
    else
        '?'
    end
end

function tounicode(p::Piece)::Char
    chars = ['♙', '♘', '♗', '♖', '♕', '♔', '?', '?', '♟', '♞', '♝', '♜', '♛', '♚']
    if isok(p)
        chars[p.val]
    else
        '?'
    end
end

"""
    isslider(t::PieceType)
    isslider(p::Piece)

Determine whether a piece is a sliding piece.
"""
function isslider(t::PieceType)::Bool
    t >= BISHOP && t <= QUEEN
end

function isslider(p::Piece)::Bool
    isslider(ptype(p))
end


text:
module ManifoldHMC

using Documenter
using SimpleTraits
using LinearAlgebra
using Markdown: @doc_str
using Manifolds:
    Manifold,
    IsDecoratorManifold,
    Euclidean,
    Sphere,
    ProductManifold,
    base_manifold,
    similar_result,
    project_point!,
    project_tangent!,
    exp!
using ForwardDiff
using DiffResults

include("flow.jl")

include("manifolds.jl")

end # module


text:
function iterate_interpolations(source::String)
    result = Union{Expr, JSString, Symbol}[]
    lastidx = 1; i = 1; lindex = lastindex(source)
    while true
        c = source[i]
        if c == '$'
            if !isempty(lastidx:(i - 1))
                push!(result, JSString(source[lastidx:(i - 1)]))
            end
            expr, i2 = Meta.parse(source, i + 1, greedy = false, raise = false)
            if i2 >= lindex && expr === nothing
                error("Invalid interpolation at index $(i)-$(lindex): $(source[i:lindex])")
            end
            i = i2
            push!(result, esc(expr))
            lastidx = i
            i > lindex && break
        else
            if i == lindex
                if !isempty(lastidx:lindex)
                    push!(result, JSString(source[lastidx:lindex]))
                end
                break
            end
            i = Base.nextind(source, i)
        end
    end
    return result
end

macro js_str(js_source)
    value_array = :([])
    append!(value_array.args, iterate_interpolations(js_source))
    return :(JSCode($value_array))
end

append_source!(x::JSCode, value::String) = push!(x.source, JSString(value))
append_source!(x::JSCode, value::JSString) = push!(x.source, value)
append_source!(x::JSCode, value::JSCode) = append!(x.source, value.source)


text:
"""
```
struct Model <: AbstractModel
    objective::Objective
    eq_constraints::VectorOfFunctions
    ineq_constraints::VectorOfFunctions
    box_min::AbstractVector
    box_max::AbstractVector
    init::AbstractVector
    integer::AbstractVector
end
```

The `Model` structs stores information about the nonlinear optimization problem.
- `objective`: the objective function of the problem of type [`Objective`](@ref).
- `eq_constraints`: the equality constraints of the problem of type [`VectorOfFunctions`](@ref). Each function in `ineq_constraints` can be an instance of [`IneqConstraint`](@ref) or [`AbstractFunction`](@ref). If the function is not an `IneqConstraint`, its right-hand-side bound is assumed to be 0.
- `ineq_constraints`: the inequality constraints of the problem of type [`VectorOfFunctions`](@ref). Each function in `ineq_constraints` can be an instance of [`IneqConstraint`](@ref) or [`AbstractFunction`](@ref). If the function is not an `IneqConstraint`, its right-hand-side bound is assumed to be 0.
"""
mutable struct Model{Tv <: AbstractVector} <: AbstractModel
    objective::Union{Nothing, Objective}
    eq_constraints::VectorOfFunctions
    ineq_constraints::VectorOfFunctions
    box_min::Tv
    box_max::Tv
    init::Tv
    integer::BitVector
end

"""
    Model()
    Model(f)

Constructs an empty model or a model with objective function `f`. The decision variables are assumed to be of type `Vector{Any}`.
"""
Model(f::Union{Nothing, Function} = nothing) = Model(Any, f)

"""
    Model(::Type{T}, f::Union{Nothing, Function}) where {T}

Constructs an empty model with objective function `f` and decision variable value type `Vector{T}`. `f` can be an instance of `Base.Function` but must return a number,  or it can be an intance of [`Objective`](@ref). `f` can also be `nothing` in which case, the objective function can be defined later using [`set_objective!`](@ref).
"""
Model(::Type{T}, f::Function) where {T} = Model(T, Objective(f))
function Model(::Type{T}, obj::Union{Nothing, Objective}) where {T}
    return Model(obj, VectorOfFunctions(EqConstraint[]), VectorOfFunctions(IneqConstraint[]), T[], T[], T[], falses(0))
end

getobjective(m::AbstractModel) = m.objective

getineqconstraints(m::AbstractModel) = m.ineq_constraints

geteqconstraints(m::AbstractModel) = m.eq_constraints

function getobjectiveconstraints(m::AbstractModel)
    return VectorOfFunctions((getobjective(m), getineqconstraints(m), geteqconstraints(m)))
end

getineqconstraint(m::AbstractModel, i::Integer) = getineqconstraints(m).fs[i]

geteqconstraint(m::AbstractModel, i::Integer) = geteqconstraints(m).fs[i]

getnineqconstraints(m::AbstractModel) = length(getineqconstraints(m))

getneqconstraints(m::AbstractModel) = length(geteqconstraints(m))

getnconstraints(m::AbstractModel) = getnineqconstraints(m) + getneqconstraints(m)

getnvars(m::AbstractModel) = length(getmin(m))

"""
    getdim(m::AbstractModel)

Returns a 2-tuple of the number of constraints and the number of variables in the model `m`.
"""
getdim(m::AbstractModel) = (getnconstraints(m), getnvars(m))

getmin(m::AbstractModel)= m.box_min
getmin(m::AbstractModel, i) = getmin(m)[i]
isinteger(m::AbstractModel, i) = m.integer[i]

getmax(m::AbstractModel) = m.box_max
getmax(m::AbstractModel, i) = getmax(m)[i]

function getinit(model::AbstractModel)
    _model, _, unflatten = tovecmodel(model)
    return unflatten(getinit(_model))
end

function setmin!(m::AbstractModel, min)
    getmin(m) .= min
    return m
end
function setmin!(m::AbstractModel, i, min)
    getmin(m)[i] = min
    return m
end

function setmax!(m::AbstractModel, max)
    getmax(m) .= max
    return m
end
function setmax!(m::AbstractModel, i, max)
    getmax(m)[i] = max
    return m
end

# Box constraints
function setbox!(m::AbstractModel, minb, maxb)
    getmin(m) .= minb
    getmax(m) .= maxb
    return m
end
function setbox!(m::AbstractModel, i, minb, maxb)
    getmin(m)[i] = minb
    getmax(m)[i] = maxb
    return m
end
function setinteger!(m::AbstractModel, i, integer)
    m.integer[i] = integer
    return m
end

function addvar!(m::Model, lb, ub; init = deepcopy(lb), integer = false)
    push!(getmin(m), lb)
    push!(getmax(m), ub)
    push!(m.init, init)
    push!(m.integer, integer)
    return m
end
function addvar!(m::Model, lb::Vector, ub::Vector; init = deepcopy(lb), integer = falses(length(lb)))
    append!(getmin(m), lb)
    append!(getmax(m), ub)
    append!(m.init, init)
    append!(m.integer, integer)
    return m
end

"""
    set_objective!(m::AbstractModel, f::Function)

Sets the objective of the moodel `m` to the function `f`. `f` must return a scalar.
"""
function set_objective!(m::AbstractModel, f::Function; kwargs...)
    m.objective = Objective(f; kwargs...)
    return m
end

function add_ineq_constraint!(m::AbstractModel, f::Function, s = 0.0; dim = length(flatten(f(getinit(m)))[1]), kwargs...)
    return add_ineq_constraint!(m, FunctionWrapper(f, dim), s; kwargs...)
end
function add_ineq_constraint!(m::AbstractModel, f::AbstractFunction, s = 0.0; kwargs...)
    return add_ineq_constraint!(m, IneqConstraint(f, s; kwargs...))
end
function add_ineq_constraint!(m::AbstractModel, f::IneqConstraint)
    push!(m.ineq_constraints.fs, f)
    return m
end
function add_ineq_constraint!(m::AbstractModel, fs::Vector{<:IneqConstraint})
    append!(m.ineq_constraints.fs, fs)
    return m
end

function add_eq_constraint!(m::AbstractModel, f::Function, s = 0.0; dim = length(flatten(f(getinit(m)))[1]), kwargs...)
    return add_eq_constraint!(m, FunctionWrapper(f, dim), s; kwargs...)
end
function add_eq_constraint!(m::AbstractModel, f::AbstractFunction, s = 0.0; kwargs...)
    return add_eq_constraint!(m, EqConstraint(f, s; kwargs...))
end
function add_eq_constraint!(m::AbstractModel, f::EqConstraint)
    push!(m.eq_constraints.fs, f)
    return m
end
function add_eq_constraint!(m::AbstractModel, fs::Vector{<:EqConstraint})
    append!(m.eq_constraints.fs, fs)
    return m
end