text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -37.92418043260588 <= q <= 19.457860872935036)

                     add_NL_constraint(m, :(softplus(-0.5948345021955626 + 0.7909419235855366*softplus(0.16572281076834727 + 0.7701313689304001*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.9350053029495258*softplus(0.6292157071886879 + 0.7447123643459479*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.6696024550400765*softplus(-0.5341808090129341 + 0.4032770914369115*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.9733717899660395*softplus(0.43846326345255404 + -0.02862602361460853*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.07261217388131413*softplus(0.02496868545483233 + -0.04485855107115411*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) + softplus(-0.7044745597894595 + -0.85619817917634*softplus(0.16572281076834727 + 0.7701313689304001*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.37830231556965854*softplus(0.6292157071886879 + 0.7447123643459479*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.12156626196118481*softplus(-0.5341808090129341 + 0.4032770914369115*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.3744591738356289*softplus(0.43846326345255404 + -0.02862602361460853*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.4985780257939494*softplus(0.02496868545483233 + -0.04485855107115411*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) + softplus(-0.9683275825894269 + -0.7637278723892034*softplus(0.16572281076834727 + 0.7701313689304001*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.9143661573707975*softplus(0.6292157071886879 + 0.7447123643459479*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.7824888723979511*softplus(-0.5341808090129341 + 0.4032770914369115*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.6866950641423624*softplus(0.43846326345255404 + -0.02862602361460853*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.09633196432882363*softplus(0.02496868545483233 + -0.04485855107115411*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) + softplus(-0.5143076265809325 + -0.11312975790997548*softplus(0.16572281076834727 + 0.7701313689304001*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.7215302240627146*softplus(0.6292157071886879 + 0.7447123643459479*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.8436266175918532*softplus(-0.5341808090129341 + 0.4032770914369115*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.8029468027318871*softplus(0.43846326345255404 + -0.02862602361460853*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.07363029133667531*softplus(0.02496868545483233 + -0.04485855107115411*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) + softplus(-0.8626839052591899 + -0.6726505975051698*softplus(0.16572281076834727 + 0.7701313689304001*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.3960198610694281*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.006655112331988011*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.5604852483433853*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.5611932928423773*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + -0.7431517048783043*softplus(0.6292157071886879 + 0.7447123643459479*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.8071073342022652*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.7192221993944883*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1911667575027427*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.30728485701247665*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.03587478002549904*softplus(-0.5341808090129341 + 0.4032770914369115*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + -0.6061337759397625*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + 0.5534397384266061*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.16208042832895408*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.7376214833764658*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.7126331133171049*softplus(0.43846326345255404 + -0.02862602361460853*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.5029397949622734*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.20686891235628035*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.9777195341190259*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.2597965376833278*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2]))) + 0.7945114361960401*softplus(0.02496868545483233 + -0.04485855107115411*softplus(0.4458559142916436 + -0.1548130095620479*$(x[1]) + 0.5123092826047779*$(x[2])) + 0.47774406442351136*softplus(0.5904352841714982 + -0.7465945594038805*$(x[1]) + -0.5086687599184327*$(x[2])) + -0.5903487696420049*softplus(-0.6551208030223536 + -0.40103431116018706*$(x[1]) + -0.6543625861790279*$(x[2])) + -0.1719483424118926*softplus(-0.8605921328376098 + -0.6449201156921633*$(x[1]) + -0.9058095618930082*$(x[2])) + -0.8023455576856944*softplus(-0.4569139457142626 + -0.5316223845746624*$(x[1]) + 0.2651376452711478*$(x[2])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
// 2 instead of 1 arguments

int main() {
	int x = foo(1,2);
	return 0 ;
}

int foo(int y) {
 return y;
}

ext: jl

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
# LP programming using MATLAB-like construction
#
#   maximize 1000 x + 350 y
#
#       s.t. x >= 30, y >= 0
#            x - 1.5y >= 0  (i.e. -x + 1.5 y <= 0)
#            12 x + 8 y <= 1000
#            1000 x + 300 y <= 70000
#
#   solution: (59.0909, 36.3636)
#   objv: 71818.1818
#

using MathProgBase
using Xpress, Compat.Test
using Compat.SparseArrays

@testset "Basics 2" begin
    model = xpress_model(
        name="lp_02",
        sense=:maximize,
        f = [1000., 350.],
        A = [-1. 1.5; 12. 8.; 1000. 300.],
        b = [0., 1000., 70000.],
        lb = [0., 30.])

    @test Xpress.get_obj(model) == [1000, 350]
    @test Xpress.get_constrmatrix(model) == sparse([-1 1.5; 12 8; 1000 300])
    @test Xpress.get_rhs(model) == [0, 1000, 70000]
    @test Xpress.get_lb(model) == [0, 30]
    @test Xpress.get_ub(model) == [1e20, 1e20]
    @test Xpress.get_sense(model) == [:<=,:<=,:<=]

    optimize(model)

    ans = get_optiminfo(model)
    @test ans.status_lp == :optimal

    @test isapprox(get_solution(model), [59.0909 , 36.3636]; atol = 1e-3)

    @test isapprox(get_objval(model), 71818.1818181818; atol = 1e-3)
end


text:
using Test
using TypeClasses
using DataTypesBasic
using Suppressor

splitln(str) = split(strip(str), "\n")

# Combine
# =======

a = State(s -> ("hello $s", s+1))
b = State(s -> ("!", s+10))

@test (a ⊕ b)(4) == ("hello 4!", 15)

# FunctorApplicativeMonad
# =======================

g = State(s -> (s*2, s+1))
f = State(s -> (s*s, s))

# just function composition
fg = map(f, g)
@test fg(3) == ((36, 6), 4)

fPLUSg = mapn(f, g) do x1, x2
  x1 + x2
end

@test fPLUSg(3) == (9+6, 4)

fPRODg = @syntax_flatmap begin
  x1 = f
  x2 = g
  @pure x1 * x2
end

@test fPRODg(3) == (9*6, 4)

putget = @syntax_flatmap begin
  putstate(4)
  x = getstate
  @pure x * x
end
@test putget(()) == (16, 4)

# FlipTypes
# =========

# this only works because of the general `combine` implementation for State
# so be cautious as the eltypes need to support `combine` to not get a MethodError in runtime

a = State.([s -> (s, s), s -> (2s, s*s), s -> (s*s, 4)])
a2 = [map(x -> pure(Vector, x), v) for v in a]
@test flip_types(a)(3) == ([3, 6, 81], 4)


text:
using JuMP
using GLPK
model = Model(GLPK.Optimizer)
#unregister(model, :x)
@variable(model,1 <= x <= 100,Int)
@variable(model,1 <= y <= 100,Int)

@NLobjective(model,Min,abs(y/x-2.9))
print(model)
optimize!(model)
solution_summary(model)
println("** Optimal objective function value = ", JuMP.objective_value(model))
println("** Optimal x = ", JuMP.value(x))
println("** Optimal y = ", JuMP.value.(y))

max_stars_repo_name: powerzf/VehicleCalculate

max_issues_repo_name: powerzf/VehicleCalculate

max_forks_repo_name: powerzf/VehicleCalculate


text:
using Distributions, DataFrames, DataFramesMeta, Arrow, Statistics, Dates, CSV, ODBC, WRDSMerger
using Revise
using AbnormalReturns

##
conn = ODBC.Connection("wrds-pgdata-64")
##

df_ff_data = @chain ff_data(conn) begin
    @rtransform(:mkt = :mktrf + :rf)
end

ff_data_max = maximum(df_ff_data.date)

@time df_crsp_raw = @chain crsp_data(conn, Date(1990), ff_data_max; cols=["ret"]) begin
    sort([:permno, :date])
end

data = MarketData(
    df_ff_data,
    df_crsp_raw
)

##

mkt_mat = @chain df_ff_data begin
    select([:mktrf, :smb, :hml, :umd])
    dropmissing
    Matrix
    cov
end
# 4×4 Matrix{Float64}:
#   0.000116248  -1.03659e-5   1.1843e-5   -1.04678e-5
#  -1.03659e-5    3.49562e-5  -2.37018e-6   1.9103e-6
#   1.1843e-5    -2.37018e-6   3.77871e-5  -1.05367e-5
#  -1.04678e-5    1.9103e-6   -1.05367e-5   5.98879e-5

##

@time df_all_regs = @chain df_crsp_raw begin
    select([:permno, :date])
    dropmissing
    @by(
        :permno,
        :date_start = minimum(:date),
        :date_end = maximum(:date)
    )
    @transform(:reg = AbnormalReturns.vector_reg(data, :permno, :date_start, :date_end, @formula(ret ~ mktrf + smb + hml + umd); minobs=100, save_residuals=true))
end

##

function winsorize(vect::AbstractArray, low_per::AbstractFloat=0.01, high_per::AbstractFloat=0.99)
    lower = quantile(skipmissing(vect), low_per)
    upper = quantile(skipmissing(vect), high_per)
    returnVect = vect[:]
    for (i, v) in enumerate(vect)
        if typeof(v) == Missing
            continue
        end
        if v > upper
            returnVect[i] = upper
        elseif v < lower
            returnVect[i] = lower
        end
    end
    return returnVect
end

df_coefs = @chain df_all_regs begin
    @rtransform(:r2 = r2(:reg))
    dropmissing
    @rtransform(:coefs = NamedTuple{(Symbol.(coefnames(:reg))..., ^(:var))}((coef(:reg)..., var(:reg))))
     _.coefs
    DataFrame
    rename("(Intercept)" => "int")
    @transform(
        :int = winsorize(:int),
        :mktrf = winsorize(:mktrf),
        :smb = winsorize(:smb),
        :hml = winsorize(:hml),
        :umd = winsorize(:umd),
        :var = winsorize(:var)
    )
end
##
coef_means = mean(df_coefs[:, 1:5] |> Matrix, dims=1)[1, :]
# [
#     0.0003416652538979692,
#     0.7557301603393168,
#     0.6035365140401592,
#     0.13841355976862427,
#     -0.0917780093857371,
# ]
coef_cov = cov(df_coefs[:, 1:5] |> Matrix)
# 2.52326e-6  -4.71181e-5   3.31182e-5  -7.05013e-5    1.32415e-5
# -4.71181e-5   0.27552      0.173033     0.0742338    -0.0448335
#  3.31182e-5   0.173033     0.389159     0.0460528    -0.0280572
# -7.05013e-5   0.0742338    0.0460528    0.40226      -0.000262376
#  1.32415e-5  -0.0448335   -0.0280572   -0.000262376   0.183439

fit_mle(Gamma, @rsubset(df_coefs, :var > 0).var)
# Gamma{Float64}(α=0.5792772418365709, θ=0.00448646291687775)


text:
using BinaryBuilder
using Pkg

name = "Darknet_CUDA"

include("../common.jl")

version, sources, script, products, dependencies = gen_common(; gpu = true)

push!(sources,
      [
          ArchiveSource("https://github.com/JuliaGPU/CUDABuilder/releases/download/v0.3.0/CUDNN+CUDA10.1.v7.6.5.x86_64-linux-gnu.tar.gz", "79de5b5085a33bc144b87028e998a1d295a15c3424d6d45b25defe500f616974", unpack_target = "cudnn"),
          DirectorySource("./bundled"),
      ]...)

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc, compiler_abi = CompilerABI(cxxstring_abi = :cxx11)),
]

push!(dependencies, BuildDependency(PackageSpec(name="CUDA_full_jll")))

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using Rasters, Adapt, Test
using Rasters.GeoFormatTypes, Rasters.LookupArrays

struct CustomArray{T,N} <: AbstractArray{T,N}
    arr::Array
end

CustomArray(x::Array{T,N}) where {T,N} = CustomArray{T,N}(x)
Adapt.adapt_storage(::Type{<:CustomArray}, xs::Array) = CustomArray(xs)

Base.size(x::CustomArray, y...) = size(x.arr, y...)
Base.getindex(x::CustomArray, y...) = getindex(x.arr, y...)

@testset "adapt AbstractProjected" begin
    l = Projected([1:10...];
        crs=ProjString("+proj="), mappedcrs=EPSG(4326),
        span=Regular(1), metadata=Metadata(:a=>"1", :b=>"2")
    )
    l1 = Adapt.adapt(CustomArray, l)
    @test parent(parent(l1)) isa CustomArray
    @test parent(parent(l1)).arr == [1:10...]
    @test metadata(l1) == NoMetadata()
    l = Mapped([1:10...];
        crs=ProjString("+proj="), mappedcrs=EPSG(4326),
        span=Regular(1), metadata=Metadata(:a=>"1", :b=>"2")
    )
    l1 = Adapt.adapt(CustomArray, l)
    @test parent(parent(l1)) isa CustomArray
    @test parent(parent(l1)).arr == [1:10...]
    @test metadata(l1) == NoMetadata()
end


text:
const SomeOf{AType} = NTuple{N,AType} where {N}

function flatten_unrolled()
    ()
end
function flatten_unrolled(container, rest...)
    container..., flatten_unrolled(rest...)...
end

function my_flatten(containers)
    # TODO: avoid recursion for large tuples
    flatten_unrolled(containers...)
end

@inline function all_unrolled()
    true
end
@inline function all_unrolled(item1, rest...)
    if item1
        all_unrolled(rest...)
    else
        false
    end
end

function my_all(items)
    # TODO: avoid recursion for large tuples
    all_unrolled(items...)
end

function map_unrolled(call, variables::Tuple{})
    ()
end
function map_unrolled(call, variables)
    call(first(variables)), map_unrolled(call, tail(variables))...
end
function map_unrolled(call, variables1::Tuple{}, variables2::Tuple{})
    ()
end
function map_unrolled(call, variables1, variables2)
    call(first(variables1), first(variables2)),
    map_unrolled(call, tail(variables1), tail(variables2))...
end
function map_unrolled(call, variables1, ::Tuple{})
    error("Mismatch in map: $variables1 left over")
end
function map_unrolled(call, ::Tuple{}, variables2)
    error("Mismatch in map: $variables2 left over")
end

function my_map(call, variables)
    # TODO: avoid recursion for large tuples
    map_unrolled(call, variables)
end
function my_map(call, variables1, variables2)
    # TODO: avoid recursion for large tuples
    map_unrolled(call, variables1, variables2)
end

function partial_map_unrolled(call, fixed, variables::Tuple{})
    ()
end
function partial_map_unrolled(call, fixed, variables)
    call(fixed, first(variables)), partial_map_unrolled(call, fixed, tail(variables))...
end
function partial_map_unrolled(call, fixed, variables1::Tuple{}, variables2::Tuple{})
    ()
end
function partial_map_unrolled(call, fixed, variables1, variables2)
    call(fixed, first(variables1), first(variables2)),
    partial_map_unrolled(call, fixed, tail(variables1), tail(variables2))...
end
function partial_map_unrolled(call, fixed, variables1, ::Tuple{})
    error("Mismatch in partial map: $variables1 left over")
end
function partial_map_unrolled(call, fixed, ::Tuple{}, variables2)
    error("Mismatch in partial map: $variables2 left over")
end

function partial_map(call, fixed, variables)
    # TODO: avoid recursion for large tuples
    partial_map_unrolled(call, fixed, variables)
end
function partial_map(call, fixed, variables1, variables2)
    # TODO: avoid recursion for large tuples
    partial_map_unrolled(call, fixed, variables1, variables2)
end

@inline in_unrolled(needle) = false
@inline in_unrolled(needle, hay, stack...) =
    if needle === hay
        true
    else
        in_unrolled(needle, stack...)
    end

setdiff_unrolled(less) = ()
function setdiff_unrolled(less, first_more, mores...)
    rest = setdiff_unrolled(less, mores...)
    if in_unrolled(first_more, less...)
        rest
    else
        first_more, rest...
    end
end

function my_setdiff(more, less)
    # TODO: avoid recursion for large tuples
    setdiff_unrolled(less, more...)
end

"""
    key(pair)

The `key` in a `key => value` `pair`.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> @inferred key(:a => 1)
:a

julia> @inferred key((:a, 1))
:a
```
"""
function key((a_key, a_value))
    a_key
end
function key(pair::Pair)
    pair.first
end
export key

"""
    value(pair)

The `value` in a `key => value` `pair`.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> @inferred value(:a => 1)
1

julia> @inferred value((:a, 1))
1
```
"""
function value((a_key, a_value))
    a_value
end
function value(pair::Pair)
    pair.second
end
export value

"""
    @if_known(something)

If `something` is `missing`, return `missing`, otherwise, `something`.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> function test(x)
           first(@if_known(x))
       end;

julia> @inferred test((1, 2))
1

julia> @inferred test(missing)
missing
```
"""
macro if_known(something)
    quote
        something = $(esc(something))
        if something === missing
            return missing
        else
            something
        end
    end
end
export @if_known


text:
using Dates

export days, yearfraction

#export AbstractDayCountConvention

#export DayCount30_360, DayCountAct,  DayCountAct_Act

export DayCount30A_360, DayCount30E_360, DayCount30E_360ISDA
export DayCountAct_360, DayCountAct_364, DayCountAct_365Fixed, DayCountAct_365Actual, DayCountAct_365Long, DayCountAct_365_25
export DayCountAct_ActISDA

abstract type AbstractDayCountConvention end

abstract type DayCount30_360 <: AbstractDayCountConvention end
abstract type DayCountAct <: AbstractDayCountConvention end
abstract type DayCountAct_Act <: AbstractDayCountConvention end

abstract type DayCount30A_360 <: DayCount30_360 end
abstract type DayCount30E_360 <: DayCount30_360 end
abstract type DayCount30E_360ISDA <: DayCount30_360 end

abstract type DayCountAct_360 <: DayCountAct end
abstract type DayCountAct_364 <: DayCountAct end
abstract type DayCountAct_365Fixed <: DayCountAct end
abstract type DayCountAct_365_25 <: DayCountAct end
abstract type DayCountAct_365Actual <: DayCountAct end
abstract type DayCountAct_365Long <: DayCountAct end

abstract type DayCountAct_ActISDA <: DayCountAct_Act end

@inline daysFormula(y1::Int, y2::Int, m1::Int, m2::Int, d1::Int, d2::Int)::Int = 360(y2-y1) + 30(m2-m1) + (d2-d1)
"""
    function days(startDate::Date, endDate::Date, ::Type{T}) where T <: DayCount30A_360

Calculates days between startDate and endDate using 2006 ISDA definitions 4.16f.
also known as 30/360 U.S. Municipal, 30/360 Bond Basis, 30/360 ISDA
"""
days(startDate::Date, endDate::Date, ::Type{T}) where T <: DayCount30A_360 = begin
    d1 = min(day(startDate), 30)
    de = day(endDate)
    d2 = if d1 == 30
        min(de, 30)
    else
        de
    end
    daysFormula(year(startDate), year(endDate), month(startDate), month(endDate), d1, d2)
end
days(startDate::Date, endDate::Date, ::Type{T}) where T <: DayCount30E_360 = begin
    daysFormula(
        year(startDate), year(endDate),
        month(startDate), month(endDate),
        min(day(startDate), 30), min(day(endDate), 30))
end
days(startDate::Date, endDate::Date, terminationDate::Date, ::Type{T}) where T <: DayCount30E_360ISDA = begin
    d1 = if startDate == Dates.lastdayofmonth(startDate)
        30
    else
        day(startDate)
    end

    d2 = if endDate == Dates.lastdayofmonth(endDate) && !(month(endDate) == 2 && endDate == terminationDate)
        30
    else
        day(endDate)
    end
    daysFormula(year(startDate), year(endDate), month(startDate), month(endDate), d1, d2)
end
yearfraction(startDate::Date, endDate::Date, ::Type{T}) where T <: DayCount30_360 =
    days(startDate, endDate, T) / 360.0
yearfraction(startDate::Date, endDate::Date, terminationDate::Date, ::Type{T}) where T <: DayCount30_360 = days(startDate, endDate, terminationDate, T) / 360.0

days(startDate::Date, endDate::Date, ::Type{T}) where T <: DayCountAct =
    (endDate - startDate).value
yearfraction(startDate::Date, endDate::Date, ::Type{T}) where T <: DayCountAct_360 =
    days(startDate, endDate, T) / 360.0
yearfraction(startDate::Date, endDate::Date, ::Type{T}) where T <: DayCountAct_364 =
    days(startDate, endDate, T) / 364.0
yearfraction(startDate::Date, endDate::Date, ::Type{T}) where T <: DayCountAct_365Fixed =
    days(startDate, endDate, T) / 365.0
yearfraction(startDate::Date, endDate::Date, ::Type{T}) where T <: DayCountAct_365Actual = begin
    d = days(startDate, endDate, T)
    for y ∈ year(startDate):year(endDate)
        if Dates.isleapyear(y) && startDate < Date(y,2,29) ≤ endDate
            return d / 366.0
        end
    end
    d / 365.0
end

max_stars_repo_name: binnisb/JuliaLearn

max_issues_repo_name: binnisb/JuliaLearn

max_forks_repo_name: binnisb/JuliaLearn


text:
# Base64 Decoder
# ==============

struct Base64Decoder <: Codec
    table::CodeTable64
    state::State
    buffer::Buffer
end

Base64Decoder(table::CodeTable64) = Base64Decoder(table, State(), Buffer(3))

"""
    Base64Decoder(;urlsafe::Bool=false, ignore::String="$(escape_string(whitespace))")

Create a base64 decoding codec.

Arguments
---------
- `urlsafe`: use `-` and `_` as the last two values
- `ignore`: ASCII characters that will be ignored while decoding
"""
function Base64Decoder(;urlsafe::Bool=false, ignore::String=whitespace)
    if urlsafe
        table = BASE64_URLSAFE
    else
        table = BASE64_STD
    end
    table = copy(table)
    ignorechars!(table, ignore)
    return Base64Decoder(table)
end

const Base64DecoderStream{S} = TranscodingStream{Base64Decoder,S} where S<:IO

"""
    Base64DecoderStream(stream::IO; kwargs...)

Create a base64 decoding stream (see `Base64Decoder` for `kwargs`).
"""
function Base64DecoderStream(stream::IO; kwargs...)
    return TranscodingStream(Base64Decoder(;kwargs...), stream)
end

function TranscodingStreams.startproc(
        codec :: Base64Decoder,
        state :: Symbol,
        error :: Error)
    start!(codec.state)
    return :ok
end

function TranscodingStreams.process(
        codec  :: Base64Decoder,
        input  :: Memory,
        output :: Memory,
        error  :: Error)
    table = codec.table
    state = codec.state
    buffer = codec.buffer

    # Check if we can encode data.
    if !is_running(state)
        error[] = ArgumentError("decoding is already finished")
        return 0, 0, :error
    elseif output.size < 3
        # Need more output space.
        return 0, 0, :ok
    end

    # Load the first bytes.
    i = j = 0
    while buffer.size < 3 && i < input.size
        buffer[buffer.size+=1] = input[i+=1]
    end
    c1 = c2 = c3 = c4 = BASE64_CODEIGN
    if buffer.size ≥ 1
        c1 = decode(table, buffer[1])
    end
    if buffer.size ≥ 2
        c2 = decode(table, buffer[2])
    end
    if buffer.size ≥ 3
        c3 = decode(table, buffer[3])
    end
    empty!(buffer)

    # Start decoding loop.
    status = :ok
    @inbounds while true
        if c1 > 0x3f || c2 > 0x3f || c3 > 0x3f || c4 > 0x3f
            i, j, status = decode64_irregular(table, c1, c2, c3, c4, input, i, output, j, error)
        else
            output[j+1] = c1 << 2 | c2 >> 4
            output[j+2] = c2 << 4 | c3 >> 2
            output[j+3] = c3 << 6 | c4
            j += 3
        end
        if i + 4 ≤ input.size && j + 3 ≤ output.size && status == :ok
            c1 = decode(table, input[i+1])
            c2 = decode(table, input[i+2])
            c3 = decode(table, input[i+3])
            c4 = decode(table, input[i+4])
            i += 4
        else
            break
        end
    end

    # Epilogue.
    if status == :end || status == :error
        finish!(state)
    end
    return i, j, status
end

# Decode irregular code (e.g. non-alphabet, padding, etc.).
function decode64_irregular(table, c1, c2, c3, c4, input, i, output, j, error)
    # Skip ignored chars.
    while true
        if c1 == BASE64_CODEIGN
            c1, c2, c3 = c2, c3, c4
        elseif c2 == BASE64_CODEIGN
            c2, c3 = c3, c4
        elseif c3 == BASE64_CODEIGN
            c3 = c4
        elseif c4 == BASE64_CODEIGN
            # pass
        else
            break
        end
        if i + 1 ≤ input.size
            c4 = decode(table, input[i+=1])
        else
            c4 = BASE64_CODEEND
            break
        end
    end

    # Write output.
    if c1 ≤ 0x3f && c2 ≤ 0x3f && c3 ≤ 0x3f && c4 ≤ 0x3f
        output[j+=1] = c1 << 2 | c2 >> 4
        output[j+=1] = c2 << 4 | c3 >> 2
        output[j+=1] = c3 << 6 | c4
        status = :ok
    elseif c1 ≤ 0x3f && c2 ≤ 0x3f && c3 ≤ 0x3f && c4 == BASE64_CODEPAD
        c4 = 0x00
        output[j+=1] = c1 << 2 | c2 >> 4
        output[j+=1] = c2 << 4 | c3 >> 2
        status = :end
    elseif c1 ≤ 0x3f && c2 ≤ 0x3f && c3 == c4 == BASE64_CODEPAD
        c3 = c4 = 0x00
        output[j+=1] = c1 << 2 | c2 >> 4
        status = :end
    elseif c1 == c2 == c3 == BASE64_CODEIGN && c4 == BASE64_CODEEND
        status = :end
    else
        error[] = DecodeError("invalid base64 data")
        status = :error
    end
    return i, j, status
end


text: ncap2 -O -v -s '*dz=dz_cT(:, 0, 0); *area=area_sT(0, :, :); *QT_vol[$time,$z_t,$nlat,$nlon]=QFLX_TEMP*dz*area; *vol[$z_t,$nlat,$nlon]=dz*area; sum_vol=vol.ttl(); QT_ttl=QT_vol.ttl($z_t, $nlat, $nlon); QT_ttll=QT_ttl.ttl(); mean_QT=QT_ttll/sum_vol;' frc.nc int_frc.nc

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
@inline function _foldl_lazy_cat_vectors(rf, acc, vectors)
    isempty(vectors) && return complete(rf, acc)
    result = @return_if_reduced foldlargs(acc, vectors...) do acc, arr
        foldl_nocomplete(rf, acc, arr)
    end
    return complete(rf, result)
end

"""
    _foldl_lazy_hcat(rf, acc, coll::LazyArrays.Hcat)
"""
@inline _foldl_lazy_hcat(rf, acc, coll::AbstractMatrix) =
    _foldl_lazy_cat_vectors(rf, acc, coll.args)
# Hcat currently always is an `AbstractMatrix`

"""
    _foldl_lazy_vcat(rf, acc, coll::LazyArrays.Vcat)
"""
@inline function _foldl_lazy_vcat(rf, acc, coll)
    isempty(coll.args) && return complete(rf, acc)
    coll isa AbstractVector && return _foldl_lazy_cat_vectors(rf, acc, coll.args)
    coll :: AbstractMatrix
    for j in axes(coll, 2)
        vectors = view.(coll.args, Ref(:), j)
        acc = @return_if_reduced _foldl_lazy_cat_vectors(rf, acc, vectors)
    end
    return complete(rf, acc)
end
# Vcat currently always is an `AbstractVector` or `AbstractMatrix`

# TODO: write reduce for Vcat/Hcat which can be done in the "natural" order


text:
# code borrow from https://github.com/JuliaPackaging/Scratch.jl/blob/v1.1.0/test/utils.jl
maybe_io(io) = (; io = io)

function temp_pkg_dir(fn::Function; rm=true)
    old_load_path = copy(LOAD_PATH)
    old_depot_path = copy(DEPOT_PATH)
    old_home_project = Base.HOME_PROJECT[]
    old_active_project = Base.ACTIVE_PROJECT[]
    try
        empty!(LOAD_PATH)
        empty!(DEPOT_PATH)
        Base.HOME_PROJECT[] = nothing
        Base.ACTIVE_PROJECT[] = nothing
        withenv("JULIA_PROJECT" => nothing,
                "JULIA_LOAD_PATH" => nothing,
                "JULIA_PKG_DEVDIR" => nothing) do
            env_dir = mktempdir()
            depot_dir = mktempdir()
            try
                push!(LOAD_PATH, "@", "@v#.#", "@stdlib")
                push!(DEPOT_PATH, depot_dir)
                Pkg.develop(PackageSpec(path=dirname(@__DIR__)); maybe_io(pkgio)...)
                fn(env_dir)
            finally
                try
                    rm && Base.rm(env_dir; force=true, recursive=true)
                    rm && Base.rm(depot_dir; force=true, recursive=true)
                catch err
                    # Avoid raising an exception here as it will mask the original exception
                    println(Base.stderr, "Exception in finally: $(sprint(showerror, err))")
                end
            end
        end
    finally
        empty!(LOAD_PATH)
        empty!(DEPOT_PATH)
        append!(LOAD_PATH, old_load_path)
        append!(DEPOT_PATH, old_depot_path)
        Base.HOME_PROJECT[] = old_home_project
        Base.ACTIVE_PROJECT[] = old_active_project
    end
end

function with_active_project(fn, project)
    old_project = Base.ACTIVE_PROJECT[]
    Base.ACTIVE_PROJECT[] = abspath(project)
    try
        fn()
    finally
        Base.ACTIVE_PROJECT[] = old_project
    end
end

function temp_project_file(uuid = nothing)
    project = joinpath(mktempdir(), "Project.toml")
    touch(project)
    if uuid !== nothing
        open(project, "w") do io
            println(io, "uuid = \"$(uuid)\"")
        end
    end
    return project
end


text:
using MusicXML
using Documenter

# using Pkg
# pkg"activate .."
# push!(LOAD_PATH,"../src/")

makedocs(;
    modules=[MusicXML],
    authors="Amin Yahyaabadi",
    repo="https://github.com/JuliaMusic/MusicXML.jl/blob/{commit}{path}#L{line}",
    sitename="MusicXML.jl",
    format=Documenter.HTML(;
        canonical="https://JuliaMusic.github.io/MusicXML.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaMusic/MusicXML.jl",
)


text:
using Bokeh
using Dates # this is for julia < 0.4, for >=0.4 you would need Base.Dates

start = now()
days = 720
x = map(d -> start + Dates.Day(d), 1:days)
y1 = cumsum(rand(days) - 0.5)
small_bit = (maximum(y1) - minimum(y1))*0.02
y2 = y1 + (2*rand(days) - 1) * small_bit
glyphs = [Glyph(:Line, linecolor="red"),
		  Glyph(:Circle, fillcolor="rgba(255, 77, 255, 0.4)", size=6, linecolor="transparent")]
plot(x, [y1 y2], glyphs, title="Share Price Prediction", legends=["Average", "Closing Price"])
showplot()


text:
# TODO Re-design test cases
# TODO In the paper they used 3000s of simulations for tests
# TODO Firing rates in the paper: 0.76Hz, 1.27Hz, 2.5Hz
# TODO For bursts, N=8
# TODO delta t could be se to 50 ms 9according to Demas et al

Δ1 = 20 # Millisecond
Δ2 = 200 # Millisecond
Δ3 = 500 # Millisecond

inter1 = 8.# Millisecond
inter2 = 10.# Millisecond

time_resolution = 1000

simplest_train = vcat([ [inter1 * k for k in 1:5] .+ Δ1,
                        [inter2 * k for k in 1:3] .+ Δ2,
                        ]...
                    ) ./time_resolution
simplest_train2 = simplest_train .+ 0.003#Second

simplest_train_doubled = vcat([
                                ([inter1/2 * k for k in 1:5] .+ Δ1 .+12),
                                ([inter2/2 * k for k in 1:3] .+ Δ2 .+10),
                                ]...
                            ) ./time_resolution
simplest_train_doubled2  = simplest_train_doubled  .+ 0.003#Second

spike_vector = vcat([ [inter1 * k for k in 1:5] .+ Δ1,
                        [inter1 * k for k in 1:3] .+ Δ2,
                        [inter2 * k for k in 1:5] .+ Δ3,
                        ]...
                    ) ./time_resolution
spike_vector2 = spike_vector .+ 0.003
Δt = 4/time_resolution
start_time = 0
stop_time = 500/ time_resolution
stop_time = 1000/ time_resolution
tailing_coefficient(spike_vector, spike_vector, Δt, start_time, stop_time,)

# Necessary properties, described in the paper introducing the measure
# N1 Symmetry: The measure C should be symmetric: for spike trains
# A and B, C(A, B)==C(B, A)
@test tailing_coefficient(spike_vector, spike_vector2, Δt, start_time, stop_time,) ==
            tailing_coefficient(spike_vector2, spike_vector, Δt, start_time, stop_time,)

# N2 Robust to variations in the firing rate: For instance, given two spike
# trains with a particular correlational structure, if the rates of both trains
# are doubled but the structure is preserved, the correlation measure should
# remain the same.

@test tailing_coefficient(
            simplest_train,
            simplest_train2,
            Δt,
            start_time,
            stop_time,
        ) == tailing_coefficient(
            simplest_train_doubled,
            simplest_train_doubled2,
            Δt,
            start_time,
            stop_time,
        ) broken=true

# N3 Robust to amount of data: In practice, this often means robust to recording
# duration.
@test abs(tailing_coefficient(
                simplest_train,
                simplest_train2,
                Δt,
                start_time,
                stop_time,
            ) - tailing_coefficient(
                simplest_train,
                simplest_train2,
                Δt,
                start_time,
                stop_time + 0.05,
            ),
        ) < 0.01

# N4 Bounded: The measure should be bounded taking a value of1 when the spike
# trains are identical, with a value of zero corresponding to no correlation and
# 1 corresponding to anti-correlation.
@test tailing_coefficient(
                simplest_train,
                simplest_train,
                Δt,
                start_time,
                stop_time,
            )  == 1 broken=true

@test tailing_coefficient(
                simplest_train,
                simplest_train.+0.5,
                Δt,
                start_time,
                stop_time,
            )  == 0 broken=true
# N5 Robust to variations in t: Small variations to t should not introduce
# artefacts into the measure.
@test abs(tailing_coefficient(
                simplest_train,
                simplest_train2,
                Δt,
                start_time,
                stop_time,
            ) - tailing_coefficient(
                simplest_train,
                simplest_train2,
                Δt+0.0005,
                start_time,
                stop_time,
            ),
        ) < 0.00108

# N6 Anticorrelation: The measure should discriminate between no correlation and
# anticorrelation.
# TODO this one I have no idea how to approach


text:
module KomePackage

function greet(arg = "Test")
    return "Hello $arg"
end

end # module


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

mutable struct UpdateResource <: SwaggerModel
    tags::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: tags

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

const _property_map_UpdateResource = Dict{Symbol,Symbol}(Symbol("tags")=>Symbol("tags"))
const _property_types_UpdateResource = Dict{Symbol,String}(Symbol("tags")=>"Dict{String, String}")
Base.propertynames(::Type{ UpdateResource }) = collect(keys(_property_map_UpdateResource))
Swagger.property_type(::Type{ UpdateResource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_UpdateResource[name]))}
Swagger.field_name(::Type{ UpdateResource }, property_name::Symbol) =  _property_map_UpdateResource[property_name]

function check_required(o::UpdateResource)
    true
end

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


text:
using GSA_CME
using Test

@testset "GSA_CME.jl" begin
    # Write your tests here.
end


text:
abstract type AbstractScalingFunction end

struct InteriorScalingFunction <: AbstractScalingFunction
    values::OffsetArrays.OffsetVector{Float64, Vector{Float64}}
    support::Vector{DyadicRational}
    filter::InteriorFilter
    resolution::Int64
end

function Base.show(io::IO, phi::InteriorScalingFunction)
    print(io, "Interior scaling function with $(vanishing_moments(phi)) vanishing moments at resolution $(resolution(phi))")
end

Base.values(phi::AbstractScalingFunction) = phi.values
filter(phi::InteriorScalingFunction) = phi.filter
vanishing_moments(phi::InteriorScalingFunction) = phi |> filter |> vanishing_moments
resolution(phi::AbstractScalingFunction) = phi.resolution
support_boundaries(phi::InteriorScalingFunction) = support_boundaries(filter(phi))
support(phi::AbstractScalingFunction) = phi.support

function find_index(phi::AbstractScalingFunction, x::DyadicRational)
    x_scale = resolution(x)
    phi_scale = resolution(phi)

    if (phi_scale < x_scale)
        throw(DomainError(x, "Scaling function not computed at this resolution"))
    end

    idx = numerator(x) << (phi_scale - x_scale)
end

function get_value(phi::AbstractScalingFunction, x::DyadicRational)
    idx = find_index(phi, x)

    v = values(phi)
    if checkbounds(Bool, v, idx)
        return v[idx]
    else
        return 0.0
    end
end

function set_value!(phi::AbstractScalingFunction, x::DyadicRational, value::Float64)
    idx = find_index(phi, x)

    if checkbounds(Bool, values(phi), idx)
        values(phi)[idx] = value
    else
        throw(BoundsError())
    end
end

function set_value!(phi::AbstractScalingFunction, x::Integer, value::Float64)
    set_value!(phi, DyadicRational(x, 0), value)
end

function get_value(phi::AbstractScalingFunction, x::Integer)
    get_value(phi, DyadicRational(x, 0))
end

function (phi::InteriorScalingFunction)(x::DyadicRational)
    get_value(phi, x)
end

function (phi::InteriorScalingFunction)(x::DyadicRational, J::Integer, k::Integer = 0)
end

function initialize_interior_scaling_function(f::InteriorFilter, R::Integer)
    support_left, support_right = support_boundaries(f)
    support = all_dyadic_rationals(support_left, support_right, R)

    phi2 = InteriorScalingFunction(
        OffsetArrays.OffsetVector{Float64}(undef, support_left*2^R:support_right*2^R),
        support, f, R
    )
end

"""
    interior_scaling_function(::InteriorFilter) -> InteriorScalingFunction

Compute an interior scaling function in the integers in its support.
"""
function interior_scaling_function(h::InteriorFilter)
    if vanishing_moments(h) == 1
        # Haar
        return InteriorScalingFunction(
            OffsetArrays.OffsetVector([1.0; 0.0], 0:1),
            DyadicRational.(0:1, 0), h, 0
        )
    end

    H = dyadic_dilation_matrix(h)

    nonzero_function_values = eigval1(H)
    # Normalize scaling function in L2
    LinearAlgebra.rmul!(nonzero_function_values, 1/sum(nonzero_function_values))
    function_values = [0.0 ; nonzero_function_values ; 0.0]

    InteriorScalingFunction(
        OffsetArrays.OffsetVector(function_values, support(h)),
        DyadicRational.(support(h), 0), h, 0
    )
end

"""
    interior_scaling_function(::InteriorFilter, ::Integer) -> InteriorScalingFunction

Compute an interior scaling function in the dyadic rationals of scale `R` in its support.
"""
function interior_scaling_function(h::InteriorFilter, R::Integer)
    if R < 0
        throw(DomainError(R, "Resolution must be non-negative"))
    end

    phi = interior_scaling_function(h)

    for _ = 1:R
        phi = increase_resolution(phi)
    end

    return phi
end

"""
    increase_resolution(::InteriorScalingFunction)

Increase the resolution of a DaubScaling scaling function by one.
"""
function increase_resolution(phi::InteriorScalingFunction)
    h = filter(phi)
    phi2 = initialize_interior_scaling_function(h, resolution(phi) + 1)
    support2 = support(phi2)

    for (index, x) in enumerate(support2)
        if isodd(index)
            set_value!(phi2, x, phi(x))
            continue
        end

        phi_val = 0.0
        for j in support(h)
            phi_val += h[j] * phi(2*x - j)
        end

        set_value!(phi2, x, sqrt2 * phi_val)
    end

    return phi2
end

"""
    dyadic_dilation_matrix(::InteriorFilter) -> Matrix

Return the "dyadic dilation matrix" from an interior filter used to compute the associated scaling
function in the integers of its support.
"""
function dyadic_dilation_matrix(h::InteriorFilter)
    sz = length(h) - 2
    left_support = support_boundaries(h)[1]

    [sqrt2 * h[2*i - j + left_support] for i = 1:sz, j = 1:sz]
end

function Base.collect(phi::AbstractScalingFunction)
    x = phi |> support .|> float
    y = phi |> values |> parent

    return x, y
end

@recipe function f(phi::AbstractScalingFunction)
    seriestype := :path

    collect(phi)
end


text:
module Day02

using AdventOfCode2019

function day02(input::String = readInput(joinpath(@__DIR__, "..", "data", "day02.txt")))
    data = parse.(Int, split(input, ","))
    dataCopy = copy(data)

    data[1+1] = 12
    data[2+1] = 2
    part1 = _run_program(data)

    copy!(data, dataCopy)
    part2 = 0
    for noun = 0:99
        for verb = 0:99
            copy!(data, dataCopy)
            data[2] = noun
            data[3] = verb
            result = _run_program(data)
            if result == 19690720
                part2 = 100 * noun + verb
            end
        end
    end
    return [part1, part2]
end

function _run_program(data::Array{Int, 1})
    for i = 1:4:length(data)
        optcode = data[i]
        if optcode == 99
            break
        elseif optcode == 1 || optcode == 2
            pos1 = data[i+1] + 1
            pos2 = data[i+2] + 1
            pos3 = data[i+3] + 1
            if optcode == 1
                data[pos3] = data[pos1] + data[pos2]
            elseif optcode == 2
                data[pos3] = data[pos1] * data[pos2]
            end
        else
            println("Invalid optcode")
            break
        end
    end
    return data[1]
end

end # module


text:
function auth_resp_amqplain(auth_params::Dict{String,Any})
    params = Dict{String,Any}("LOGIN" => auth_params["LOGIN"], "PASSWORD" => auth_params["PASSWORD"])
    iob = IOBuffer()
    write(iob, TAMQPFieldTable(params))
    bytes = take!(iob)
    skipbytes = sizeof(fieldtype(TAMQPFieldTable, :len))
    bytes = bytes[(skipbytes+1):end]
    TAMQPLongStr(bytes)
end

const AUTH_PROVIDERS = Dict{String,Function}("AMQPLAIN" => auth_resp_amqplain)


text:
# Drive train example
# ============================

using HybridSystems, MathematicalSystems, LazySets, Plots
import SX:is_halfspace, is_hyperplane, convert, parse_sxmath

x1 = parse_sxmath("x1 - -0.0488 == 0.0005600000000000001 * (x4 - 20.0)")
x2 = parse_sxmath("x2 - -15.67 == 0.46699999999999997 * (x4 - 20.0)")
x3 = parse_sxmath("x3 == 0.0")
x4 = parse_sxmath("x4 <= 40.0")
x5 = parse_sxmath("x5 == 0.0")
x6 = parse_sxmath("x6 - 20.0 == 1.0 * (x4 - 20.0)")
x7 = parse_sxmath("x7 - 240.0 == 12.0 * (x4 - 20.0)")
x8 = parse_sxmath("x8 - -0.0019199999999999998 == 0.00005999999999999998 * (x4 - 20.0)")
x9 = parse_sxmath("x9 - 20.0 == 1.0 * (x4 - 20.0)")
t = parse_sxmath("t == 0.0")

expr_p = parse_sxmath("x1 - -0.0488 == 0.0005600000000000001 * (x4 - 20.0) & x2 - -15.67 == 0.46699999999999997 * (x4 - 20.0) & x3 == 0.0 & x5 == 0.0 & x6 - 20.0 == 1.0 * (x4 - 20.0) & x7 - 240.0 == 12.0 * (x4 - 20.0) & x8 - -0.0019199999999999998 == 0.00005999999999999998 * (x4 - 20.0) & x9 - 20.0 == 1.0 * (x4 - 20.0) & 20.0 <= x4 & x4 <= 40.0 & t == 0.0")
expr = [x1, x2, x3, x4,x5,x6,x7,x8,x9,t]
x0sets = []
vars = Basic[:x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :t]
for expr_i in expr_p
    if is_halfspace(expr_i)
        push!(x0sets, convert(HalfSpace, expr_i, vars=vars))
    elseif is_hyperplane(expr_i)
        push!(x0sets, convert(Hyperplane, expr_i, vars=vars))
    end
end

AFFINE_SYSTEM = ConstrainedLinearControlContinuousSystem
HS = readsxmodel("models/ARCH2018/Drivetrain/drivetrain_theta1_5percent_flat_manually.xml", raw_dict=true)
# initial condition in mode 1

# calculate reachable states up to time T
prob = InitialValueProblem(HS, x0sets)
input_options = Options(:mode=>"reach")

problem_options = Options(:vars=>[1,3,10], :T=>10.0, :δ=>0.005, :plot_vars=>[1, 3], :verbosity=>1);
options_input = merge(problem_options, input_options)
using Polyhedra
sol = solve(HS, options_input);

plot(sol, indices=1:2:length(sol.Xk))

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
using OpenCV
using ArgParse
using Cambrian
using CartesianGeneticProgramming
using IICGP
import Cambrian.mutate  # mutate function scope definition

function load_img(rom_name::String, frame_number::Int64)
    filename = string(@__DIR__, "/images/", rom_name, "_frame_$frame_number.png")
    return OpenCV.imread(filename)
end

function generate_io_image(rom_name::String="freeway", frame_number::Int64=30)
    img = load_img(rom_name, frame_number)
    r, g, b = IICGP.split_rgb(img)

    # Arbitrary application of simple OpenCV functions
    i = IICGP.CGPFunctions.f_add_img(r, g)
    j = IICGP.CGPFunctions.f_erode_img(i, i)
    k = IICGP.CGPFunctions.f_compare_eq_img(j, g)
    l = IICGP.CGPFunctions.f_dilate_img(k, k)
    target = IICGP.CGPFunctions.f_compare_ge_img(j, l)

    return [r, g, b], target
end

function fitness(ind::CGPInd, input::Vector{T}, target::T) where {T <: OpenCV.InputArray}
    output = process(ind, input)
    [-OpenCV.norm(output[1], target)]
end

# Read configuration file
s = ArgParseSettings()
@add_arg_table! s begin
    "--cfg"
    help = "configuration script"
    default = "cfg/ipcgp.yaml"
    "--seed"
    help = "random seed for evolution"
    arg_type = Int
    default = 0
    "--ind"
    help = "individual for evaluation"
    arg_type = String
    default = ""
end
args = parse_args(ARGS, s)
n_in = 3  # RGB images
n_out = 1  # Single image
inp, target = generate_io_image()
img_size = size(target)
cfg = CartesianGeneticProgramming.get_config(args["cfg"]; function_module=IICGP.CGPFunctions, n_in=n_in, n_out=n_out, img_size=img_size)

#=
# Test I/O without evolution
foo = IPCGPInd(cfg)
process(foo, inp)
# IICGP.display_buffer(foo, 2, indexes=1:3)
=#

if length(args["ind"]) > 0
    ind = CGPInd(cfg, read(args["ind"], String))
    ftn = fitness(ind, inps, outs)
    println("Fitness: ", ftn)
else
    # Define mutate and fit functions
    mutate(i::CGPInd) = goldman_mutate(cfg, i, init_function=IPCGPInd)
    fit(i::CGPInd) = fitness(i, inp, target)
    # Create an evolution framework
    e = CGPEvolution(cfg, fit, init_function=IPCGPInd)
    # Run evolution
    run!(e)
end


text:
dotestmean   = true
doprintmeans = false
reltol       = .01          # required test accuracy
algs      = (Direct(),SortingDirect())

# run the given number of SSAs and return the mean
function runSSAs(jump_prob, Nsims, idx)
    Psamp = zeros(Int, Nsims)
    for i in 1:Nsims
        sol = solve(jump_prob, SSAStepper())
        Psamp[i] = sol[idx,end]
    end
    mean(Psamp)
end

function execute_test(u0, tf, rates, rs, Nsims, expected_avg, idx, test_name)
    prob = DiscreteProblem(u0, (0.0, tf), rates)

    for method in algs
        jump_prob = JumpProblem(prob, method, rs, save_positions=(false,false))
        avg_val = runSSAs(jump_prob, Nsims, idx)

        if dotestmean
            if doprintmeans
                println(test_name, ", method = ", typeof(method), ", mean = ", avg_val, ", act_mean = ", expected_avg)
            end
            @test abs(avg_val - expected_avg) < reltol * expected_avg
        end
    end
end

# nonlinear reaction test
Nsims = 32000
tf = .01
u0 = [200, 100, 150]
expected_avg = 84.876015624999994
rs = @reaction_network dtype begin
    k1, 2A --> B
    k2, B --> 2A
    k3, A + B --> C
    k4, C --> A + B
    k5, 3C --> 3A
end k1 k2 k3 k4 k5
rates = [1., 2., .5, .75, .25]
execute_test(u0, tf, rates, rs, Nsims, expected_avg, 1, "Nonlinear rx test")

# DNA repression model
rs = @reaction_network ptype begin
    k1, DNA --> mRNA + DNA
    k2, mRNA --> mRNA + P
    k3, mRNA --> 0
    k4, P --> 0
    k5, DNA + P --> DNAR
    k6, DNAR --> DNA + P
end k1 k2 k3 k4 k5 k6
Nsims        = 8000
tf           = 1000.0
u0           = [1,0,0,0]
expected_avg = 5.926553750000000e+02
rates = [.5, (20*log(2.)/120.), (log(2.)/120.), (log(2.)/600.), .025, 1.]
execute_test(u0, tf, rates, rs, Nsims, expected_avg, 3, "DNA test")

# DNA repression model, mix of jump types
rs = @reaction_network ptype begin
    k1*DNA, 0 --> mRNA
    k2*mRNA, 0 --> P
    k3, mRNA --> 0
    k4, P --> 0
    k5, DNA + P --> DNAR
    k6, DNAR --> DNA + P
end k1 k2 k3 k4 k5 k6
Nsims        = 8000
tf           = 1000.0
u0           = [0,0,0,0]
u0[something(findfirst(isequal(:DNA),rs.syms),0)] = 1
expected_avg = 5.926553750000000e+02
rates = [.5, (20*log(2.)/120.), (log(2.)/120.), (log(2.)/600.), .025, 1.]
execute_test(u0, tf, rates, rs, Nsims, expected_avg, something(findfirst(isequal(:P),rs.syms), 0), "DNA mixed jump test")

# simple constant production with degratation
rs = @reaction_network pdtype begin
    1000., 0 --> A
    10, A --> 0
end
rates = nothing
Nsims        = 16000
tf           = 1.0
u0           = [0]
expected_avg = 1000.0 /10*(1. - exp(-10*tf))
execute_test(u0, tf, rates, rs, Nsims, expected_avg, 1, "Zero order test")

# this is just a test to make sure a mixture of jumps actually runs without crashes
network = @reaction_network rnType  begin
    0.01, (X,Y,Z) --> 0
    hill(X,3.,100.,-4), 0 --> Y
    hill(Y,3.,100.,-4), 0 --> Z
    hill(Z,4.5,100.,-4), 0 --> X
    hill(X,2.,100.,6), 0 --> R
    hill(Y,15.,100.,4)*0.002, R --> 0
    20, 0 --> S
    R*0.005, S --> SP
    0.01, SP + SP --> SP2
    0.05, SP2 --> 0
end;
prob = DiscreteProblem([200.,60.,120.,100.,50.,50.,50.], (0.,4000.))
for method in algs
    jump_prob = JumpProblem(prob, method, network)
    sol = solve(jump_prob,SSAStepper());
end

# make sure problem instantiation works when mass action jumps
# rate consts are a nontrivial expression on the params
network = @reaction_network begin
    1, X --> 2*X
    1/K, 2X --> X
end K
p = [1000.]
prob = DiscreteProblem([500.], (0.,100.), p)
jump_prob = JumpProblem(prob, Direct(), network)

# same as above
network = @reaction_network begin
    p1*p2, X --> Y
end p1 p2
p = [1.,2.]
prob = DiscreteProblem([10.,10.],(0.,10.), p)
jump_prob = JumpProblem(prob, Direct(), network)


text:
daipath = expanduser("../../misc/samcnet/deps/libdai")
if !isfile("libdaiwrap.so") || stat("libdaiwrap.so").mtime < stat("dai.cpp").mtime
	run(`g++ dai.cpp -I$daipath/include -L$daipath/lib -ldai -lgmpxx -lm -shared -fPIC -Wl,-rpath,\$ORIGIN -o libdaiwrap.so`)
end
!isfile("libdai.so") && run(`ln -s $daipath/lib/libdai.so .`)
#push!(DL_LOAD_PATH, pwd())

#using BinDeps

#@BinDeps.setup

##libdai_base = library_dependency("libdai_base",aliases=["libdai"])
#daipath = expanduser("~/GSP/research/samc/samcnet/deps/libdai")
##provides(Binaries,joinpath(daipath,"lib"),libdai_base)

#libdai = library_dependency("libdaiwrap")

#prefix=joinpath(BinDeps.depsdir(libdai),"usr")
#provides(BuildProcess,
	#(@build_steps begin
		#FileRule(joinpath(prefix,"lib","libdaiwrap.so"),@build_steps begin
		#`g++ dai.cpp -I$daipath/include -L$daipath/lib -ldai -lgmpxx -lm -shared -fPIC -o libdaiwrap.so`
		#`mkdir -p $prefix/lib`
		#`mv libdaiwrap.so $prefix/lib`
		#`ln -s $daipath/lib/libdai.so $prefix/lib/`
		#end)
	#end),libdai)

##provides(Sources,URI("http://staff.science.uva.nl/~jmooij1/libDAI/libDAI-0.3.1.tar.gz"), libdai)
## will need to modify the libdai config file to enable JunctionTree algorithm

#@BinDeps.install


text:
export StackFrames, ResizeImage

using ImageTransformations: imresize!
import CircularArrayBuffers
using CircularArrayBuffers: CircularArrayBuffer
using MacroTools: @forward

"""
    ResizeImage(img::Array{T, N})
    ResizeImage(dims::Int...) -> ResizeImage(Float32, dims...)
    ResizeImage(T::Type{<:Number}, dims::Int...)

Using BSpline method to resize the `state` field of an observation to size of `img` (or `dims`).
"""
struct ResizeImage{T,N}
    img::Array{T,N}
end

ResizeImage(dims::Int...) = ResizeImage(Float32, dims...)
ResizeImage(T::Type{<:Number}, dims::Int...) = ResizeImage(Array{T}(undef, dims))

function (p::ResizeImage)(state::AbstractArray)
    imresize!(p.img, state)
    p.img
end

"""
    StackFrames(::Type{T}=Float32, d::Int...)

Use a pre-initialized [`CircularArrayBuffer`](@ref) to store the latest several states specified by `d`. Before processing any observation, the buffer is filled with `zero{T}
by default.
"""
struct StackFrames{T,N} <: AbstractArray{T,N}
    buffer::CircularArrayBuffer{T,N}
end

@forward StackFrames.buffer Base.size, Base.getindex
Base.IndexStyle(x::StackFrames) = IndexStyle(x.buffer)

StackFrames(d::Int...) = StackFrames(Float32, d...)

function StackFrames(::Type{T}, d::Vararg{Int,N}) where {T,N}
    p = StackFrames(CircularArrayBuffer{T}(d...))
    for _ in 1:CircularArrayBuffers.capacity(p.buffer)
        push!(p.buffer, zeros(T, size(p.buffer)[1:N-1]))
    end
    p
end

function (p::StackFrames{T,N})(state::AbstractArray) where {T,N}
    push!(p.buffer, state)
    p
end

function RLBase.reset!(p::StackFrames{T,N}) where {T,N}
    fill!(p.buffer, zero(T))
    p
end

"""
When pushing a `StackFrames` into a `CircularArrayBuffer` of the same dimension,
only the latest frame is pushed. If the `StackFrames` is one dimension lower,
then it is treated as a general `AbstractArray` and is pushed in as a frame.
"""
function Base.push!(cb::CircularArrayBuffer{T,N}, p::StackFrames{T,N}) where {T,N}
    push!(cb, select_last_frame(p.buffer))
end


text:
# Test Interpolation
@testset "lagrange order 1" begin
    order = 1
    ip = Lagrange{1,RefLine,order}()
    bf = getnbasefunctions(ip)
    nv = CGMethod1D.nvertexdofs(ip)
    @test isapprox(bf, order+1, atol=1e-8)
    @test isapprox(nv, 1, atol=1e-8)
end

@testset "lagrange order 2" begin
    order = 2
    ip = Lagrange{1,RefLine,order}()
    bf = getnbasefunctions(ip)
    nv = CGMethod1D.nvertexdofs(ip)
    @test isapprox(bf, order+1, atol=1e-8)
    @test isapprox(nv, 1, atol=1e-8)
end

@testset "value order 1" begin
    order = 1
    ip = Lagrange{1,RefLine,order}()
    i1 = 1
    i2 = 2
    ξ = Vec{1}((-1.0,))
    N1 = value(ip, i1, ξ)
    N2 = value(ip, i2, ξ)
    println("N1:$N1; N2:$N2")
    @test isapprox(N1+N2, 1, atol=1e-8)
end

@testset "value order 2" begin
    order = 2
    ip = Lagrange{1,RefLine,order}()
    i1 = 1
    i2 = 2
    i3 = 3
    ξ = Vec{1}((-1.0,))
    N1 = value(ip, i1, ξ)
    N2 = value(ip, i2, ξ)
    N3 = value(ip, i3, ξ)
    println("N1:$N1; N2:$N2; N3:$N3")
    @test isapprox(N1+N2+N3, 1, atol=1e-8)
end

@testset "reference_coordinates order 1" begin
    order = 1; n = order+1;
    ip = Lagrange{1,RefLine,order}()
    coords = reference_coordinates(ip)
    println(coords)
    ξ_vec = collect(range(-1, stop=1, length=n))
    for (idx,val) in enumerate(ξ_vec)
        @test coords[idx][1] in ξ_vec
    end
    @test coords[1][1] == -1.0
    @test coords[2][1] == 1.0
    @test issorted(coords[3:end-1])
end

@testset "reference_coordinates order 2" begin
    order = 20; n = order+1;
    ip = Lagrange{1,RefLine,order}()
    coords = reference_coordinates(ip)
    println(coords)
    ξ_vec = collect(range(-1, stop=1, length=n))
    for (idx,val) in enumerate(ξ_vec)
        @test coords[idx][1] in ξ_vec
    end
    @test coords[1][1] == -1.0
    @test coords[2][1] == 1.0
    @test issorted(coords[3:end-1])
end

@testset "Ferrite 1D order 1" begin
    order = 1
    i1 = 1
    i2 = 2
    ξ = Vec{1}((0.25,))
    ip_F = Lagrange{1,RefCube,order}()
    bf_F = getnbasefunctions(ip_F)
    ip_C = Lagrange{1,RefLine,order}()
    bf_C = getnbasefunctions(ip_C)
    N1_F = value(ip_F, i1, ξ)
    N2_F = value(ip_F, i2, ξ)
    N1_C = value(ip_C, i1, ξ)
    N2_C = value(ip_C, i2, ξ)
    @test isapprox(bf_F, bf_C, atol=1e-8)
    @test isapprox(N1_F, N1_C, atol=1e-8)
    @test isapprox(N2_F, N2_C, atol=1e-8)
end

@testset "Ferrite 1D order 2" begin
    order = 2
    i1 = 1
    i2 = 2
    i3 = 3
    ξ = Vec{1}((0.25,))
    ip_F = Lagrange{1,RefCube,order}()
    bf_F = getnbasefunctions(ip_F)
    ip_C = Lagrange{1,RefLine,order}()
    bf_C = getnbasefunctions(ip_C)
    N1_F = value(ip_F, i1, ξ)
    N2_F = value(ip_F, i2, ξ)
    N3_F = value(ip_F, i3, ξ)
    N1_C = value(ip_C, i1, ξ)
    N2_C = value(ip_C, i2, ξ)
    N3_C = value(ip_C, i3, ξ)
    @test isapprox(bf_F, bf_C, atol=1e-8)
    @test isapprox(N1_F, N1_C, atol=1e-8)
    @test isapprox(N2_F, N2_C, atol=1e-8)
    @test isapprox(N3_F, N3_C, atol=1e-8)
end


text:
# Autogenerated wrapper script for Clingo_jll for x86_64-linux-musl-cxx11
export clasp, clingo, gringo, libclingo, lpconvert, reify

JLLWrappers.@generate_wrapper_header("Clingo")
JLLWrappers.@declare_executable_product(clasp)
JLLWrappers.@declare_executable_product(clingo)
JLLWrappers.@declare_executable_product(gringo)
JLLWrappers.@declare_library_product(libclingo, "libclingo.so.4")
JLLWrappers.@declare_executable_product(lpconvert)
JLLWrappers.@declare_executable_product(reify)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        clasp,
        "bin/clasp",
    )

    JLLWrappers.@init_executable_product(
        clingo,
        "bin/clingo",
    )

    JLLWrappers.@init_executable_product(
        gringo,
        "bin/gringo",
    )

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

    JLLWrappers.@init_executable_product(
        lpconvert,
        "bin/lpconvert",
    )

    JLLWrappers.@init_executable_product(
        reify,
        "bin/reify",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# ---------------------------------------------------------------------------- #
#
#   genSTL.jl
#
#   Write STL file from lattice frame representation
#
#   intrico
#   Spring 2018
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    genSTL( mesh::MeshF,  lat::Lattice, flname::String;
            n = 8, boundflat = Vector{Int64}( 0 ) )

Generates a binary .stl file for the lattice in `mesh` with the areas defined in `lat`.
`boundflat` lists which boundaries are flat.
"""
function genSTL( mesh::MeshF,  lat::Lattice, flname::String;
                 n::Int64 = 8, boundflat = Vector{Int64}( 0 ) )

    edgWrite = Vector{Vector{Bool}}( size(mesh.e,1) )
    for jj in 1:size(mesh.e,1)
        edgWrite[jj] = [true]
    end

    # open STL file
    fid = open( flname, "w" )

    # write header
    for jj in 1:80
        write(fid,' ')
    end
    write(fid,UInt32(0)) # number of faces is unknown now, so write 0

    # generate all faces (and write)
    factot = [ 0 ]
    actptsNods, edgesNods = genFacesNodsCyl( mesh, lat, n, factot, fid, edgWrite, boundflat )

    # close STL
    close( fid )

    ## reopen to write the number of faces
    fid = open( flname, "a+" )
    seekstart(fid) # go back to top of file

    # rewrite header
    for jj in 1:80
        write(fid,' ')
    end
    # write correct number of faces
    write( fid, UInt32( factot[1] ) )
    # close STL
    close( fid )

    return actptsNods, edgesNods
end

"""
    genSTL( mesh::MeshF,  lat::Lattice, flnames::Vector{String}, edgWrite::Vector{Vector{Bool}};
            n = 8, boundflat = Vector{Int64}( 0 ) )

Generates a binary .stl file for the lattice in `mesh` with the areas defined in `lat`.
`edgWrite` lists which struts to write in which file.
`boundflat` lists which boundaries are flat.
"""
function genSTL( mesh::MeshF,  lat::Lattice, flnames::Vector{String},
                 edgWrite::Vector{Vector{Bool}}; n::Int64 = 8, boundflat = Vector{Int64}( 0 ) )

    # open STL file
    fid = Vector{IOStream}( length(flnames) )
    for ii in 1:length(flnames)

        fid[ii] = open( flnames[ii], "w" )

        # write header
        for jj in 1:80
            write(fid[ii],' ')
        end
        write(fid[ii],UInt32(0)) # number of faces is unknown now, so write 0

    end

    # generate all faces (and write)
    factot = fill( 0, length(flnames) )
    actptsNods, edgesNods = genFacesNodsCyl( mesh, lat, n, factot, fid, edgWrite, boundflat )

    # close STL
    for ii in 1:length(flnames)
        close( fid[ii] )
    end

    ## reopen to write the number of faces
    for ii in 1:length(flnames)

        fid = open( flnames[ii], "a+" )
        seekstart(fid) # go back to top of file

        # rewrite header
        for jj in 1:80
            write(fid,' ')
        end
        # write correct number of faces
        write( fid,UInt32( factot[ii] ) )
        # close STL
        close( fid )

    end

    return actptsNods, edgesNods

end

"""
    genFacesNods( mesh::MeshF, lat::Lattice, fdist::Matrix{Float64}, n::Int64, fid::IOStream )

Writes the STL facets for all nodes for ASCII files.
"""
function genFacesNodsCyl( mesh::MeshF, lat::Lattice, n::Int64, nfac::Vector{Int64},
                          fid::Union{IOStream,Vector{IOStream}},
                          edgWrite::Vector{Vector{Bool}}, boundflat::Vector{Int64} )

    # setup points to cylinder connectivity
    ptsCylEdge = Vector{Vector{Vector{SVector{3,Float64}}}}( size(mesh.e,1) )
    for kk in 1:size(mesh.e,1)
        ptsCylEdge[kk] = Vector{Vector{SVector{3,Float64}}}(2)
    end

    # setup distance vector
    fdist = Vector{Vector{Float64}}( size(mesh.e,1) )
    for jj in 1:size(mesh.e,1)
        fdist[jj] = fill( 0.0, 2 )
    end

    # get flat vector, if exists
    nvecflat = Vector{SVector{3,Float64}}( length(boundflat) )
    for jj in 1:length(boundflat)
        indf = find(mesh.f[:,5] .== -boundflat[jj])[1]
        vec1 = SVector{3,Float64}( mesh.p[ mesh.f[indf,2],: ] - mesh.p[ mesh.f[indf,1],: ] )
        vec2 = SVector{3,Float64}( mesh.p[ mesh.f[indf,3],: ] - mesh.p[ mesh.f[indf,1],: ] )
        nvecflat[jj] = cross(vec1,vec2) / norm( cross(vec1,vec2) )
    end

    # generate nodes
    actptsNods = Vector{Vector{SVector{3,Float64}}}( mesh.n )
    edgesNods  = Vector{Vector{SVector{5,Int64}}}( mesh.n )
    for nn in 1:mesh.n

        genFacesNod( nn, mesh, lat, n, ptsCylEdge, nfac, fid,
                     edgWrite, fdist, boundflat, nvecflat, actptsNods, edgesNods )

    end

    # Collision warnings
    collisionWarning( mesh, fdist )

    # generate cylinders
    for kk in 1:size(mesh.e,1)
        if lat.ar[ kk ] > 0.0 && sqrt( lat.ar[ kk ] / π ) > 1.e-14
            genFacesCyl( kk, mesh, ptsCylEdge[kk], nfac, fid, edgWrite[kk] )
        end
    end

    return actptsNods, edgesNods

end

"""
    compDist( mesh::MeshF,  lat::Lattice )

Computes the distance from the end face to the node
"""
function genFacesNod( kk::Int64, mesh::MeshF,  lat::Lattice, nn::Int64,
                      ptsCylEdge::Vector{Vector{Vector{SVector{3,Float64}}}},
                      nfac::Vector{Int64},
                      fid::Union{IOStream,Vector{IOStream}},
                      edgWrite::Vector{Vector{Bool}},
                      fdist::Vector{Vector{Float64}},
                      boundflat::Vector{Int64}, nvecflat::Vector{SVector{3,Float64}},
                      actptsNods::Vector{Vector{SVector{3,Float64}}},
                      edgesNods::Vector{Vector{SVector{5,Int64}}} )

    # current coordinate
    pcurr = SVector{3}( mesh.p[kk,:] )

    # initialize connectivity information
    actptsNods[kk] = Vector{SVector{3,Float64}}( 0 )
    edgesNods[kk]  = Vector{SVector{5,Int64}}( 0 )

    # generate pairs of edges
    edg   = [i for i in 1:length(mesh.n2e[kk])]
    #   check for min area
    inddel = Vector{Int64}(0)
    for ii in 1:length(mesh.n2e[kk])
        if lat.ar[ mesh.n2e[kk][ii] ] < 0.0 || sqrt( lat.ar[ mesh.n2e[kk][ii] ] / π ) < 1.e-14
            append!(inddel,ii)
        end
    end
    edg0     = [i for i in 1:length(edg)-length(inddel)]
    indedges = [jj for jj in 1:length(edg)]
    deleteat!(indedges,inddel)

    if length(edg0) == 0 # if all nodes are zero thickness, can skip the whole function
        return nothing
    end

    edg1  =   edg0'
    eedg  =   edg0 .+ 0*edg1
    eedg1 = 0*edg0 .+   edg1
    pairs = hcat( eedg[:], eedg1[:] )

    pairs = unique(sort(pairs,2), 1) # unique pairs
    pairs = pairs[ sortperm( pairs[:,1] + pairs[:,2]*100 ), : ] # ensure highest indices are last
    # TODO: don't need to do this for each loop, can just precompute this for n = 40 and then just take part of the unique pairs
    rmax = sqrt( maximum( lat.ar[ mesh.n2e[kk] ] ) / π )

    # compute crossing point and normal vector
    normvec  = Vector{Vector{SVector{3,Float64}}}( length(edg) )
    eveclink = Vector{Vector{Int64}}( length(edg) )
    # tangvec  = Vector{Vector{Vector{Float64}}}( length(edg) )
    lvmax = fill( 0.0, length(edg) )
    for jj in 1:length(edg)
        normvec[jj]  = Vector{SVector{3,Float64}}( length(edg0)-1 )
        eveclink[jj] = Vector{Int64}( length(edg0)-1 )
        # tangvec[jj]  = Vector{Vector{Float64}}( length(edg)-1 )
    end
    indcnt = fill(1, length(edg) )
    for jj in 1:size(pairs,1)
        edg_i1 = indedges[ pairs[jj,1] ]
        edg_i2 = indedges[ pairs[jj,2] ]
        if edg_i1 == edg_i2
            continue
        end
        e1 = mesh.n2e[kk][ edg[ edg_i1 ] ]
        e2 = mesh.n2e[kk][ edg[ edg_i2 ] ]

        # compute angle between the edges
        nod1 = mesh.e[e1, find( mesh.e[e1,:] .!= kk )[] ]
        nod2 = mesh.e[e2, find( mesh.e[e2,:] .!= kk )[] ]
        vec1 = mesh.p[ nod1, : ] - mesh.p[ kk, : ]
        vec2 = mesh.p[ nod2, : ] - mesh.p[ kk, : ]
        vec1 ./= norm(vec1)
        vec2 ./= norm(vec2)
        arg = dot(vec1,vec2)
        arg = minimum( ( 1.0, arg) ) # guard for round-off errors
        arg = maximum( (-1.0, arg) ) # guard for round-off errors
        θ = acos( arg )

        # compute absolute distance between node center and intersection
        lv = rmax / tan( θ/2 ) # offset such that rods do not collide
        l = sqrt( lv^2 + rmax^2 ) # absolute distance
        lvmax[edg_i1] = max( lvmax[edg_i1], lv )
        lvmax[edg_i2] = max( lvmax[edg_i2], lv )

        magdiff = norm(vec1 - vec2)
        magave  = norm(vec1 + vec2)

        # normal vectors
        normvec[edg_i1][ indcnt[edg_i1] ] = SVector{3}( (vec1-vec2) / magdiff )
        normvec[edg_i2][ indcnt[edg_i2] ] = SVector{3}( (vec2-vec1) / magdiff )

        # add edge links
        eveclink[edg_i1][ indcnt[edg_i1] ] = e2
        eveclink[edg_i2][ indcnt[edg_i2] ] = e1

        indcnt[edg_i1] += 1
        indcnt[edg_i2] += 1

    end

    # add offset to rods
    for ii in 1:length(lvmax)
        lengthedge = norm( mesh.p[ mesh.e[mesh.n2e[kk][edg[ii]],1],: ] - mesh.p[ mesh.e[mesh.n2e[kk][edg[ii]],2],: ] )
        if lvmax[ii] < 0.30 * lengthedge # 0.15
            lvmax[ii] += 0.10 * lengthedge # 0.05
        end

        # add to global distance vector
        ie = find( mesh.e[ mesh.n2e[kk][ edg[ii] ], 1:2] .!= kk )[]
        fdist[ mesh.n2e[kk][ edg[ii] ] ][ ie ] = lvmax[ii]
    end

    # check if on boundary that has to be flat
    notflat = true
    nflat = 0
    for eecnt in 1:length(edg0)
        ee = indedges[ eecnt ]
        e1 = mesh.n2e[kk][ edg[ ee ] ]
        if any(mesh.e[ e1, 3 ] .== boundflat) || any(mesh.e[ e1, 4 ] .== boundflat)
            # TEMP for now assume we only case about one flat surface
            # append!( normvec[ee], nvecflat[1:1] )
            notflat = false
            nflat   = 2
        end
    end

    # loop over each edge
    for eecnt in 1:length(edg0)
        ee = indedges[ eecnt ]

        e1     = mesh.n2e[kk][ edg[ ee ] ]
        nod    = mesh.e[e1, find( mesh.e[e1,:] .!= kk )[] ]
        evec   = mesh.p[ nod, : ] - mesh.p[ kk, : ]
        evec ./= norm(evec)

        nnormals = length(normvec[ee]) # or should this be eecnt?

        # find intersections between planes and check if they are active
        testpt   = mesh.p[ nod, : ] - mesh.p[ kk, : ] # distance to this point is measured to determine which is closer, distance is relative to current node
        nunique  = uniqueNumPairs( nnormals )
        intersec = Vector{SVector{3,Float64}}( max( 2*(nunique - nnormals) + 1 + nflat, 3 ) )
        # edge_pairs = Vector{SVector{2,Int64}}( max( 2*(nunique - nnormals) + 1 + nflat, 3 ) )

        nact = 0
        if nunique == 1 # if nunique is 1, there are no intersections

            indother    = find( indedges .!= ee )[]
            eother      = mesh.n2e[kk][ edg[ indedges[ indother ] ] ]
            nodother    = mesh.e[eother, find( mesh.e[eother,:] .!= kk )[] ]
            evecother   = SVector{3}( (mesh.p[ nodother, : ] - mesh.p[ kk, : ]) / norm(mesh.p[ nodother, : ] - mesh.p[ kk, : ]) )

            lvec = evec + evecother
            if norm(lvec) < 1e-14
                lvec = normvec[ee][1]
            end

            lvec2 = lvec * rmax / norm( lvec - dot( lvec, evec )*evec )

            nact += 1
            intersec[nact] = SVector{3}( 1.0 * lvec2)
            # edge_pairs[nact] = SVector{2}( sort(mesh.n2e[kk][ edg[ indedges ] ]) )

            nact += 1
            intersec[nact] = SVector{3}(-1.0 * lvec2)
            # edge_pairs[nact] = SVector{2}( sort(mesh.n2e[kk][ edg[ indedges ] ]) )

        else
            for jj in 1:nunique # note: number of planes is length(edg)-1, so can reuse part of upairs

                # TODO: need to check here for alignment with normal plane, because now that is possible
                # it is needed, because point won't return active (or divide by zero), actually it probably will?
                # no doesn't work, becuase will never be active

                if pairs[jj,1] == pairs[jj,2]
                    continue
                end
                ip1 = pairs[jj,1]
                ip2 = pairs[jj,2]

                lvec     = cross( normvec[ee][ ip1 ], normvec[ee][ ip2 ] ) # always going through the node center (because both planes go through that point)
                nrmperp  = norm( lvec - dot( lvec, evec )*evec )

                # if norm( cross( normvec[ee][ ip2 ], evec ) ) < 1e-14 # problematic edge is added last, so will always be in the second column (ip2)
                #     # normal plane splits edge down the middle, need to take extra care
                #     lvec1 = lvec * rmax / nrmperp
                #
                #     lvec2 = - lvec * rmax / nrmperp
                #
                # else

                # check first point
                lvec1 = lvec * rmax / nrmperp

                act = checkPointAct( lvec1, testpt, normvec[ee], evec )

                if act && ( notflat || dot( lvec1, nvecflat[1] ) <= 1e-14 )
                    nact += 1
                    intersec[nact] = lvec1
                end

                # check second point
                lvec2 = - lvec * rmax / nrmperp
                # check ...
                act = checkPointAct( lvec2, testpt, normvec[ee], evec )
                if act && ( notflat || dot( lvec2, nvecflat[1] ) <= 1e-14 ) # last statement only evaluated if flat
                    nact += 1
                    intersec[nact] = lvec2
                end

                # end

            end
        end

        ## check for flatspots
        flatact  = fill( false, nact )
        if !notflat
            for jj in 1:length(normvec[ee])
                lvec     = cross( normvec[ee][ jj ], nvecflat[1] ) # always going through the node center (because both planes go through that point)
                nrmperp  = norm( lvec - dot( lvec, evec )*evec )

                # check first point
                lvec1 = lvec * rmax / nrmperp

                act = checkPointAct( lvec1, testpt, normvec[ee], evec )

                if act
                    nact += 1
                    intersec[nact] = lvec1
                    append!(flatact, true)
                end

                # check second point
                lvec2 = - lvec * rmax / nrmperp
                # check ...
                act = checkPointAct( lvec2, testpt, normvec[ee], evec )
                if act # last statement only evaluated if flat
                    nact += 1
                    intersec[nact] = lvec2
                    append!(flatact, true)
                end

            end

        end

        ## order points counterclockwise
        # pick first point as ϕ = 0
        zerovec = (intersec[1] - dot( intersec[1], evec ) * evec) /
                   norm( intersec[1] - dot( intersec[1], evec ) * evec )
        ϕ = fill(0.0,nact)
        for jj in 2:nact
            currnvec = intersec[jj] - dot( intersec[jj], evec ) * evec
            ϕ[jj] = acos( max( min( dot(zerovec,currnvec) / norm(currnvec), 1.0), -1.0 ) )
            if dot( cross( zerovec, currnvec ), evec ) < 0.0
                ϕ[jj] = 2*π - ϕ[jj]
            end
        end

        # sort based on ϕ
        angleind         = sortperm( ϕ )
        ϕ[1:nact]        = ϕ[angleind]
        intersec[1:nact] = intersec[angleind]
        flatact          = flatact[angleind]
        corrind          = fill( 0, nact )
        corrind[1]       = angleind[1]

        nact = 1
        for jj in 1:length(ϕ)-1
            if (2*π - ϕ[jj+1]) < 1e-14
                flatact[1] = flatact[1] || any(flatact[jj+1:end])
            else
                if (ϕ[jj+1] - ϕ[jj]) > 1e-6
                    nact += 1
                    corrind[nact] = jj+1
                else
                    flatact[corrind[nact]] = flatact[corrind[nact]] || flatact[jj+1]
                end
            end
        end

        ϕ       = ϕ[ corrind[1:nact] ]
        flatact[1] = any(flatact[1:corrind[2]-1])
        flatact = flatact[ corrind[1:nact] ]
        intersec[1:nact] = intersec[ corrind[1:nact] ]

        # add last point to array to close the loop
        if (2*π - ϕ[end]) > 1e-3
            append!( ϕ, 2*π )
            append!( flatact, flatact[1] )
            intersec[nact+1] = intersec[1]
        else
            ϕ[end] = 2*π
            nact -= 1
        end

        # add active points to overall node points
        actptlinks = [ [i for i in 1:nact]; 1 ] # if no points have been added yet
        # if ee > 1
            jj = 1 # first point
            indnod = findActPointsNod( intersec[jj], actptsNods[kk], actptlinks, jj )
            actptlinks[end] = indnod
            for jj in 2:nact
                indnod = findActPointsNod( intersec[jj], actptsNods[kk], actptlinks, jj )
            end
        # end

        # add connectivity for entire node
        edge_pairs = Vector{SVector{5,Int64}}( nact ) # [nod1,nod2,strut1,strut2]

        ## Generate points along cuts
        perpvec = cross(evec,zerovec)
        perpvec ./= norm(perpvec)

        iedge = find( mesh.e[e1,:] .== kk )[]
        ptsCylEdge[e1][iedge] = Vector{SVector{3,Float64}}(0)
        for jj in 1:nact

            # TODO need to check here if flatact is true, if true need to stitch

            # figure out on which plane to project
            ϕhalf   = ( ϕ[jj] + ϕ[jj+1] ) / 2.
            currpt  = rmax * ( zerovec * cos(ϕhalf) + perpvec * sin(ϕhalf) )
            normind, halfintersec = findNormPlane( currpt, testpt, normvec[ee], evec )

            normcurr = normvec[ee][normind]

            # Add information about which edge links the struts
            indhalf = findActPointsNod( halfintersec, actptsNods[kk] )
            findEdgesNod( edgesNods[kk], SVector{5}( [ actptlinks[jj]; actptlinks[jj+1];
                                                       sort([e1, eveclink[ee][normind]]); indhalf ] ) )

            projectCircle = true
            if flatact[jj] && flatact[jj+1]
                for qq in 1:nnormals
                    if abs( dot( nvecflat[1], normvec[ee][qq] ) ) < 1e-14
                        projectCircle = false
                        break
                    end
                end
                normcurr = - nvecflat[1]
            end

            Δϕ = ϕ[jj+1] - ϕ[jj]
            @assert Δϕ > 1e-14
            np = max( convert( Int64, ceil( (Δϕ-1e-5)/(2*π) * nn ) + 1 ), 2 ) # ensure we have at least two points

            if !projectCircle
                np = 2
            end

            ϕcurr = linspace( ϕ[jj], ϕ[jj+1], np )

            # first part
            intpts = Vector{SVector{3,Float64}}( np )
            cylpts = Vector{SVector{3,Float64}}( np )

            if projectCircle
                for ii in 1:np
                    currpt = rmax * ( zerovec * cos(ϕcurr[ii]) + perpvec * sin(ϕcurr[ii]) ) # current point in plane perpendicular to evec
                    # find actual point near node
                    d = dot( -currpt, normcurr ) / dot( evec, normcurr ) # scalar value for which line intersects plane
                    intpts[ii] = SVector{3}( currpt + d * evec + pcurr )
                    cylpts[ii] = SVector{3}( sqrt(lat.ar[e1]/π) * ( zerovec * cos(ϕcurr[ii]) + perpvec * sin(ϕcurr[ii]) ) + pcurr + evec * lvmax[ee] )

                end
            else
                intpts[1] = intersec[jj]   + pcurr
                intpts[2] = intersec[jj+1] + pcurr
                cylpts[1] = SVector{3}( sqrt(lat.ar[e1]/π) * ( zerovec * cos(ϕcurr[1]) + perpvec * sin(ϕcurr[1]) ) + pcurr + evec * lvmax[ee] )
                cylpts[2] = SVector{3}( sqrt(lat.ar[e1]/π) * ( zerovec * cos(ϕcurr[2]) + perpvec * sin(ϕcurr[2]) ) + pcurr + evec * lvmax[ee] )
            end
            # write to file
            for ii in 1:np-1
                vert = SVector{3,SVector{3,Float64}}( cylpts[ii], intpts[ii], intpts[ii+1] )
                writeFacetSTLB( vert, nfac, fid, edgWrite[e1] )

                vert = SVector{3,SVector{3,Float64}}( intpts[ii+1], cylpts[ii+1], cylpts[ii] )
                writeFacetSTLB( vert, nfac, fid, edgWrite[e1] )

            end

            # if flatact is true, stitch the node
            if flatact[jj] && flatact[jj+1]
                # write to file
                for ii in 1:np-1
                    vert = SVector{3,SVector{3,Float64}}( intpts[ii+1], intpts[ii], pcurr )
                    writeFacetSTLB( vert, nfac, fid, edgWrite[e1] )
                end
            end

            # save cylinder points for cylinder connections
            append!( ptsCylEdge[e1][iedge], cylpts[1:end-1] )
            if jj == nact
                append!( ptsCylEdge[e1][iedge], cylpts[end:end] )
            end

        end

    end

end

"""
    compDist( mesh::MeshF,  lat::Lattice )

Computes the distance from the end face to the node
"""
function genFacesCyl( cyl::Int64, mesh::MeshF,
                      ptsCylEdge::Vector{Vector{SVector{3,Float64}}}, nfac::Vector{Int64},
                      fid::Union{IOStream,Vector{IOStream}}, edgWrite::Vector{Bool} )

    frac = 0.25

    e1 = mesh.e[cyl,1]
    e2 = mesh.e[cyl,2]

    ## first edge
    # loop over first edge

    n1 = length( ptsCylEdge[1] )
    n2 = length( ptsCylEdge[2] )

    # find which side has most points
    i1 = 1
    i2 = 2
    if n2 < n1
        i2 = 1
        i1 = 2
        n1 = (n1 + n2)
        n2 = n1 - n2
        n1 = n1 - n2
    end

    ## generate the all facets, starting with the side with the fewest points
    # find first index and generate first facet
    ind0 = findMinDistInd( ptsCylEdge, i1, i2, 1, frac )
    vert = SVector{3,SVector{3,Float64}}( ptsCylEdge[i1][1], ptsCylEdge[i1][2], ptsCylEdge[i2][ind0] )

    writeFacetSTLB( vert, nfac, fid, edgWrite )

    ind00 = ind0 # need to save this index to close loop

    # generate rest of facets (also for other ring)
    for jj in 2:n1-1

        ind1 = findMinDistInd( ptsCylEdge, i1, i2, jj, frac )

        vert = SVector{3,SVector{3,Float64}}( ptsCylEdge[i1][jj], ptsCylEdge[i1][jj+1], ptsCylEdge[i2][ind1] )
        writeFacetSTLB( vert, nfac, fid, edgWrite )

        # generate other faces
        # NOTE: due to ordering ind1 < ind0
        if ind0 > ind1
            for kk = 0:(ind0-ind1-1)
                vert = SVector{3,SVector{3,Float64}}( ptsCylEdge[i2][ ind1+kk ],
                                                      ptsCylEdge[i2][ ind1+kk+1 ],
                                                      ptsCylEdge[i1][ jj ] )
                writeFacetSTLB( vert, nfac, fid, edgWrite )
            end
        elseif ind0 < ind1
            inds   = vcat(ind1:n2-1, 1:ind0-1)
            indsp1 = vcat(ind1+1:n2, 2:ind0)
            for kk in 1:length(inds)
                vert = SVector{3,SVector{3,Float64}}( ptsCylEdge[i2][ inds[kk] ],
                                                      ptsCylEdge[i2][ indsp1[kk] ],
                                                      ptsCylEdge[i1][ jj ] )
                writeFacetSTLB( vert, nfac, fid, edgWrite )
            end
        end

        ind0 = ind1

    end
    # last patch
    ind1 = ind00
    if ind0 > ind1
        for kk = 0:(ind0-ind1-1)
            vert = SVector{3,SVector{3,Float64}}( ptsCylEdge[i2][ ind1+kk ],
                                                  ptsCylEdge[i2][ ind1+kk+1 ],
                                                  ptsCylEdge[i1][ 1 ] )
            writeFacetSTLB( vert, nfac, fid, edgWrite )
        end
    elseif ind0 < ind1
        inds   = vcat(ind1:n2-1, 1:ind0-1)
        indsp1 = vcat(ind1+1:n2, 2:ind0)
        for kk in 1:length(inds)
            vert = SVector{3,SVector{3,Float64}}( ptsCylEdge[i2][ inds[kk] ],
                                                  ptsCylEdge[i2][ indsp1[kk] ],
                                                  ptsCylEdge[i1][ 1 ] )
            writeFacetSTLB( vert, nfac, fid, edgWrite )
        end
    end

end

function checkPointAct( currpt::SVector{3,Float64}, testpt::Vector{Float64},
                        normvec::Vector{SVector{3,Float64}}, evec::Vector{Float64} )

    dist    =  norm( currpt - testpt )
    distorg = dist
    act = false
    for ii in 1:length(normvec)
        # loop over intersecting planes
        d = dot( - currpt, normvec[ii] ) / dot( evec, normvec[ii] ) # scalar value for which line intersects plane
        ipt = currpt + d * evec

        cdist = norm( ipt - testpt )

        if cdist < dist
            # this point is closer
            dist = cdist
        end
    end

    # is point active? Only if new distance is same as original distance
    if abs(dist - distorg)/abs(dist) < 1e-14
        act = true
    end

    return act

end

function findNormPlane( currpt::SVector{3,Float64}, testpt::Vector{Float64},
                        normvec::Vector{SVector{3,Float64}}, evec::Vector{Float64} )

    dist = 3*norm( currpt - testpt )
    ind  = 0

    intpt = SVector{3}( 0., 0., 0. )

    for ii in 1:length(normvec)
        # loop over intersecting planes
        d = dot( - currpt, normvec[ii] ) / dot( evec, normvec[ii] ) # scalar value for which line intersects plane
        ipt = currpt + d * evec

        cdist = norm( ipt - testpt )

        if cdist < dist
            # this point is closer
            dist  = cdist
            ind   = ii
            intpt = ipt
        end
        # end
    end

    return ind, intpt

end

function findMinDistInd( pts::Vector{Vector{SVector{3,Float64}}},
                         ind1::Int64, ind2::Int64,
                         jj::Int64, frac::Float64)

    testpt = (1-frac) * pts[ind1][jj] + frac * pts[ind1][jj+1] # check for minimum distance to middle of edge

    # find minimum distance to node on other side
    dist = 1e4
    ind  = 0 # look for first index, but after that only allow for + or - one (because should be ordered)
    for kk in 1:length( pts[ind2] )
        cdist = norm( pts[ind2][kk] - testpt )
        if cdist < dist
            dist = cdist
            ind = kk
        end
    end

    return ind
end

function uniqueNumPairs( n::Int64 )
    #   (n + k - 1) nCr (k)
    # here k = 2, so
    #   (n + 1) nCr (k)
    return binomial( n + 1, 2 )
end

function collisionWarning( mesh::MeshF, fdist::Vector{Vector{Float64}} )

    nw = 0
    for jj in 1:length(fdist)

        if fdist[jj][1] > 0.0
            ledge = norm( mesh.p[ mesh.e[jj,2], : ] - mesh.p[ mesh.e[jj,1], : ] )

            if ledge < sum( fdist[jj] )
                warn( "Collision detected at edge ", jj )
                if fdist[jj][1] > ledge || fdist[jj][2] > ledge
                    println("   protrudes!")
                end
                nw += 1
            end
        end

    end

    if nw > 0
        println( "" )
        warn( "Generated ", nw, " collisions" )
    end

end

function findActPointsNod( intersec::SVector{3,Float64},
                           act_points_nod::Vector{SVector{3,Float64}},
                           actptlinks::Vector{Int64}, ind::Int64 )
    tol = 1e-10

    indnew = 0
    for jj in 1:length(act_points_nod)
        if norm( intersec - act_points_nod[jj] ) < tol
            indnew = jj
            break
        end
    end

    # add new points
    if indnew == 0
        append!( act_points_nod, [intersec] )
        indnew = length(act_points_nod)
    end

    actptlinks[ind] = indnew

    return indnew

end

function findActPointsNod( intersec::SVector{3,Float64},
                           act_points_nod::Vector{SVector{3,Float64}} )
    tol = 1e-10

    indnew = 0
    for jj in 1:length(act_points_nod)
        if norm( intersec - act_points_nod[jj] ) < tol
            indnew = jj
            break
        end
    end

    # add new points
    if indnew == 0
        append!( act_points_nod, [intersec] )
        indnew = length(act_points_nod)
    end

    return indnew

end

function findEdgesNod( edges_nod::Vector{SVector{5,Int64}}, edge::SVector{5,Int64} )

    indnew = 0
    for jj in 1:length(edges_nod)
        if edge[3] == edges_nod[jj][3] && edge[4] == edges_nod[jj][4]
            indnew = jj
            break
        end
    end

    # add new edge...
    if indnew == 0
        append!( edges_nod, [edge] )
    end

end


text:
function _JCircle(center, radius, color, linewidth, action)
    sethue(color)
    setline(linewidth)
    circle(center, radius, action)
    return center
end

"""
    1. JCircle(center::Point, radius::Real; kwargs...)
    2. JCircle(center_x::Real, center_y::Real, radius::Real; kwargs...)
    3. JCircle(p1::Point, p2::Point; kwargs...)
        - A circle that touches `p1` and `p2`
    4. JCircle(radius::Real)
        - A circle at the origin

Draw a circle at `center` with the given `radius`

# Keywords for all
- `color` = "black"
- `linewidth` = 1
- `action::Symbol` :stroke by default can be `:fill` or other actions explained in the Luxor documentation.

Returns the center of the circle
"""
JCircle(center::Point, radius::Real; color = "black", linewidth = 1, action = :stroke) =
    (
        args...;
        center = center,
        radius = radius,
        color = color,
        linewidth = linewidth,
        action = action,
    ) -> _JCircle(center, radius, color, linewidth, action)

JCircle(center_x::Real, center_y::Real, radius::Real; kwargs...) =
    JCircle(Point(center_x, center_y), radius; kwargs...)

JCircle(p1::Point, p2::Point; kwargs...) =
    JCircle(midpoint(p1, p2), distance(p1, p2) / 2; kwargs...)

JCircle(radius::Real; kwargs...) = JCircle(O, radius; kwargs...)