text:
"""
    defaultSettings()

Defines the list of configuration variables and returns their default values
in a Dict.
"""
function defaultSettings()
    # Return the default settings. All parameters must be registered here.
    Dict(
        "biggenelength" => 200, # length of the compatibility gene's sequence (if `usebiggenes`)
        "bodytemp"=> 312.0, # body temperature in K for homeothermic species
        "borders" => "absorbing", # border behaviour: absorbing/reflecting/mainland
        #"bufferoutput" => false, # use output buffering (faster execution, more RAM usage) #TODO
        "burn-in" => 1000, # timesteps before invasion starts
        "capgrowth" => false, # strictly limit individuals' size to `repsize`
        "cellsize" => 20e6, # maximum biomass per hectare in gramm (based on Clark et al. 2001)
        "compressgenes" => true, # save gene sequences as integers to reduce memory usage
        "config" => "", # configuration file name
        "debug" => false, # write out debug statements
        "degpleiotropy" => 0.1, # How frequent are pleitropy and polygenic inheritance? 0 <= degpleiotropy < 1
        "dest" => string(Dates.today()), # output folder name
        "dispmean" => 1.0, # maximum value of randomly drawn mean dispersal distance in cells
        "dispshape" => 1.0, # maximum value of randomly drawn shape parameter for dispersal kernel. determines tail fatness/long distance dispersal
        "dispmortality" => false, #boolean toggle for switching mortality dependent on dispersal distance, depends on dispfactor(see below), only in zdisperse()
        "dispfactor"=> 2, #: when dispmean multiplied by dispfactor is larger than  evolved dispersaldistance -> individual dies with 100 % chance, when smaller, chances of survival higher
        "disturbance" => 0, # percentage of individuals killed per update per cell
        "dumpindforfasta"=> false, #should dumpinds be called simultaneously with the fasta output?
        "fasta" => "off", # record fasta data? "off", "compat", "all" (high detail output)
        "fastaoutfreq" => 10, #frequency with which to write fasta data
        "fertility" => exp(28.0), # global base reproduction rate 23.8 from Brown et al. 2004, alternatively 25.0, default 30.0
        "fixtol" => true, # fix mating tolerance globally to `tolerance`
        "global-species-pool" => 0, # size of the global species pool (invasion source)
        "globalmating" => false, # global pollen "dispersal"
        "growthrate" => exp(25.2), # global base growth/biomass production from Brown et al. 2004
        "heterozygosity" => false, # keep track of heterozygosity when studying hybridisation
        "hybridisation" => "off", #toggle hybridisation on and off to allow simultaneaous hybridisation and speciation modes "ecological" and "neutral"
        "indsize" => "seed", # initialize organisms as seed, adult or mixed
        "isolationweight" => 3.0, # additional distance to be crossed when dispersing from or to isolated patches
        "lineages" => false, # record lineage and diversity data (low detail output)
        "linkage" => "random", # gene linkage type (random/full/none)
        "logging" => false, # write output to logfile
        "maps" => "", # comma-separated list of map files
        "maxbreadth" => 5.0, # maximum niche breadth
        "maxloci" => 1, # maximum number of loci/copies per gene
        "maxprec" => 10.0, # max optimum precipitation
        "maxrepsize" => 14.0, # maximal repsize in grams calculated as exp(maxrepsize) -> 1.2 t
        "maxseedsize" => 10.0, # maximal seedsize in grams calculated as exp(maxseedsize) -> 22 kg
        "maxtemp" => 313.0, # max optimum temp in K
        "mortalitytype" => "metabolic", # metabolic/habitat/global
        "minprec" => 0.0, # min optimum precipitation
        "minrepsize" => 3.0, # minimal repsize in grams calculated as exp(minrepsize) -> 20 g
        "minseedsize" => -2.0, # minimal seedsize in grams calculated as exp(minseedsize) -> 0.14 g
        "mintemp" => 273.0, # min optimum temp in K
        "mode" => "default", # experiment type ("default", "invasion", or "zosterops")
        "mortality" => exp(22), # global base mortality from Brown et al. 2004 is 26.3, but competition and dispersal introduce add. mort.
        "mutate" => true, # mutations occur
        "mutationrate" => 3.6e10, # one mutation per generation/individual, corrected for metabolic function
        "nniches" => 2, # number of environmental niches (max. 2) #XXX currently, only 2 is sensible
        "outfreq" => 10, # output frequency
        "overfill" => 1.0, # how much to overfill grid cells beyond their capacity
        "perfecttol" => 2.0, #if mortality is habitat dependent, how tolerant should a perfectly adapted individual be? scales for specialist selection pressure
        "phylconstr" => 0.1, # phylogenetic constraint during mutation and inter-loci variation. scales trait value as sd.
        "popsize" => "metabolic", # initialisation algorithm: metabolic/bodysize/minimal/single
        "precrange" => 0, # max optimum precipitation - deprecated! (use `maxprec`)
        "propagule-pressure" => 0, # number of non-native individuals introduced per invasion event
        "quiet" => false, # don't write output to screen
        "raw" => true, # record raw data
        "scaleadaptation" => 0.0, #either 0.0 or 1.0, needed parameter for the gausscurve funciton. For phylogeny experiment set to 1.0!
        "sdprec" => 0.0, # SD of precipitation change per time step
        "sdtemp" => 0.0, # SD of temperature change per time step
        "seed" => 0, # for the RNG, seed = 0 -> random seed
        "smallgenelength" => 20, # standard gene sequence length (max. 21)
        "speciation" => "neutral", # allow lineage differentiation? off/neutral/ecological
        "species" => Dict{String,Any}[], # define trait values for each Zosterops species
        "static" => false, # whether mainland sites undergo eco-evolutionary processes (implies "mainland" borders)
        "stats" => true, # record population statistics (medium detail output)
        "tolerance" => 0.8, # sequence similarity threshold for reproduction if `fixtol` == true
        # minimal required traitnames
        "traitnames" => ["compat", # place-holder for the neutral ("compatibility") gene
                         "dispmean", # mean dispersal distance
                         "dispshape", # dispersal shape parameter
                         "numpollen", # number of offspring
                         "precopt", # precipitation/resource optimum
                         "prectol", # precipitation/resource tolerance
                         "repsize", # reproductive size
                         "seqsimilarity", # required genetic sequence similarity for mating
                         "selfing", # probability of self-fertilisation
                         "seedsize", # seed size
                         "tempopt", # temperature optimum
                         "temptol"], # temperature tolerance
        "traitsforecospec" => ["dispmean", # mean dispersal distance
                               "dispshape", # dispersal shape parameter
                               "numpollen", # number of offspring
                               "precopt", # precipitation/resource optimum
                               "prectol", # precipitation/resource tolerance
                               "repsize", # reproductive size
                               "seqsimilarity", # required genetic sequence similarity for mating
                               "selfing", # probability of self-fertilisation
                               "seedsize", # seed size
                               "tempopt", # temperature optimum
                               "temptol"], # temperature tolerance
        "usebiggenes" => true, # use a longer sequence for the compatibility gene
    )
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
using DiffEqJump, DiffEqBase, Test

@time begin
  @time @testset "Constant Rate Tests" begin include("constant_rate.jl") end
  @time @testset "Variable Rate Tests" begin include("variable_rate.jl") end
  @time @testset "Monte Carlo Tests" begin include("monte_carlo_test.jl") end
  @time @testset "Split Coupled Tests" begin include("splitcoupled.jl") end
  @time @testset "SSA Tests" begin include("ssa_tests.jl") end
  @time @testset "Tau Leaping Tests" begin include("regular_jumps.jl") end
  @time @testset "SIR Discrete Callback Test" begin include("sir_model.jl") end
  @time @testset "Linear Reaction SSA Test" begin include("linearreaction_test.jl") end
  @time @testset "Mass Action Jump Tests; Gene Expr Model" begin include("geneexpr_test.jl") end
  @time @testset "Mass Action Jump Tests; Nonlinear Rx Model" begin include("bimolerx_test.jl") end
  @time @testset "Mass Action Jump Tests; Special Cases" begin include("degenerate_rx_cases.jl") end
  @time @testset "Composition-Rejection Table Tests" begin include("table_test.jl") end
  @time @testset "Extinction test" begin include("extinction_test.jl") end
  @time @testset "Saveat Regression test" begin include("saveat_regression.jl") end
  @time @testset "Ensemble Uniqueness test" begin include("ensemble_uniqueness.jl") end
  @time @testset "Thread Safety test" begin include("thread_safety.jl") end
end

ext: jl


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -8.186501529651974 <= q <= 9.622483160840265)

                     add_NL_constraint(m, :(gelu(-0.9342689023247241 + -0.06643065035242302*gelu(-0.6967397920955185 + 0.4612305646899797*gelu(-0.7651411836742423 + 0.736339264438953*$(x[1]) + 0.06323249947984699*$(x[2]) + 0.6326865446578034*$(x[3]) + 0.8935006364653741*$(x[4])) + 0.5917923558177942*gelu(0.8959549878207418 + 0.6412415470262633*$(x[1]) + 0.6383904188579423*$(x[2]) + 0.9706157764666923*$(x[3]) + 0.8108132764672078*$(x[4]))) + 0.13022424784864484*gelu(-0.8092681495291991 + -0.6577545463201151*gelu(-0.7651411836742423 + 0.736339264438953*$(x[1]) + 0.06323249947984699*$(x[2]) + 0.6326865446578034*$(x[3]) + 0.8935006364653741*$(x[4])) + 0.816846389391447*gelu(0.8959549878207418 + 0.6412415470262633*$(x[1]) + 0.6383904188579423*$(x[2]) + 0.9706157764666923*$(x[3]) + 0.8108132764672078*$(x[4])))) + gelu(0.9502905692637253 + 0.2731932353707758*gelu(-0.6967397920955185 + 0.4612305646899797*gelu(-0.7651411836742423 + 0.736339264438953*$(x[1]) + 0.06323249947984699*$(x[2]) + 0.6326865446578034*$(x[3]) + 0.8935006364653741*$(x[4])) + 0.5917923558177942*gelu(0.8959549878207418 + 0.6412415470262633*$(x[1]) + 0.6383904188579423*$(x[2]) + 0.9706157764666923*$(x[3]) + 0.8108132764672078*$(x[4]))) + 0.7525923835832526*gelu(-0.8092681495291991 + -0.6577545463201151*gelu(-0.7651411836742423 + 0.736339264438953*$(x[1]) + 0.06323249947984699*$(x[2]) + 0.6326865446578034*$(x[3]) + 0.8935006364653741*$(x[4])) + 0.816846389391447*gelu(0.8959549878207418 + 0.6412415470262633*$(x[1]) + 0.6383904188579423*$(x[2]) + 0.9706157764666923*$(x[3]) + 0.8108132764672078*$(x[4])))) - $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:
#CONSTANTS
#K = 0.005454 - converts diameter squared in square inches to square feet.
#0.00007854 converts diamter square centimeters to square meters
const K = 0.005454154
const KMETRIC = 0.00007854

abstract type Log end

mutable struct LogSegment<:Log
small_end_diam
large_end_diam
length
#shape #?
end

abstract type Shape end

function area(diameter)
    a=K * diameter^2
    return a
end

#Table 6.1 form Kershaw et al Forest Mensuration
#Equations to compute Cubic Volume of Important Solids

# Equation 6.1
#V=AₗL where Aₗ is the area of the base(large end).
mutable struct Cylinder
length
large_end_diam
end

function volume(solid::Cylinder)
    V = area(solid.large_end_diam)*solid.length
    return V
end

#Equation 6.2
#V= 1/2(AₗH)
#Paraboloid

mutable struct Paraboloid
length
large_end_diam
end

function volume(solid::Paraboloid)
    V = 1/2(area(solid.large_end_diam))*solid.length
    V
end

#Equation 6.3
#V=1/3(AₗL)
#Cone
#PREDICT CUBIC FOOT VOLUME FOR A CONIC SEGMENT SUCH AS A STEM TIP
mutable struct Cone
length
large_end_diam
end

function volume(solid::Cone)
    V = 1/3(area(solid.large_end_diam))*solid.length
    return V
end

#Equation 6.4
#V=1/4(AₗL)
#Neiloid

mutable struct Neiloid
length
large_end_diam
end

function volume(solid::Neiloid)
    V = 1/4(area(solid.large_end_diam))*solid.length
    return V
end

#Equation 6.5- Smalian's formula
#V=L/2(Aₗ + Aₛ) where L is the length and Aₛ is the area of the upper end (or small end)
#Paraboloid frustrum

#Equation 6.6 - Huber's formula
#V=AₘL where Aₘ is the area at middle point
#Paraboloid frustrum

#Equation 6.9 - Newton's formula
#V=L/6(Aₗ + 4*Aₘ + Aₛ)
#Neiloid,Paraboloid or Conic Frustrum

mutable struct ParaboloidFrustrum
length
large_end_diam
mid_point_diam #can set to nothing ( or missing in 0.7.0+?)
small_end_diam
end

function volume(solid::ParaboloidFrustrum; huber=false, newton = false)
    if huber == true
        V = area(solid.mid_point_diam) * solid.length
    elseif newton == true
        V = (solid.length/6) * (area(solid.large_end_diam) + 4*area(solid.mid_point_diam) + area(solid.small_end_diam))
    else
        V = area(solid.large_end_diam) + area(solid.small_end_diam) * (solid.length/2)
    end
return V
end

#Equation 6.7
#V=L/3(Aₗ + sqrt(Aₗ*Aₛ) + Aₛ)

mutable struct ConeFrustrum
length
large_end_diam
mid_point_diam #can set to nothing
small_end_diam
end

function volume(solid::ConeFrustrum; newton=false)
    if newton == true
        V = (solid.length/6) * (area(solid.large_end_diam) + 4*area(solid.mid_point_diam) + area(solid.small_end_diam))
        else
            V = area(solid.large_end_diam) + area(solid.small_end_diam) * (solid.length/2)
    end
return V
end

#Equation 6.8
#V=L/4(Aₗ + cbrt(Aₗ²*Aₛ) + cbrt(Aₗ*Aₛ²) + Aₛ)

mutable struct NeiloidFrustrum
length
large_end_diam
mid_point_diam #can set to nothing
small_end_diam
end

function volume(solid::NeiloidFrustrum; newton=false)
    if newton == true
        V = (solid.length/6) * (area(solid.large_end_diam) + 4*area(solid.mid_point_diam) + area(solid.small_end_diam))
    else
        V = (solid.length/4) *
        (cbrt(area(solid.large_end_diam)^2 * area(solid.small_end_diam)) +
        cbrt((K * solid.large_end_diam^2) * area(solid.small_end_diam)^2) +
        area(solid.small_end_diam))
    end
return V
end

macro LogSegment(a,b,c,d,shape)
    return :($shape($a,$b,$c,$d))
end

#@LogSegment(24,12,nothing,8,ParaboloidFrustrum)

#Function to calculate the Scribner scale volume
#for more info see: http://oak.snr.missouri.edu/forestry_functions/scribnerbfvolume.php

## Scribner Decimal C table
using OffsetArrays
scribner_decimal_c=OffsetArray([[5 5 5 5 5 5 10 10 10 10 10 10 20 20 20 20 20];
[5 5 5 5 5 5 10 10 10 10 10 10 20 20 20 20 20];
[5 5 5 10 10 10 10 20 20 20 20 20 30 30 30 30 30];
[5 10 10 10 10 10 20 20 20 20 20 20 30 30 30 30 30];
[10 10 10 20 20 20 30 30 30 30 30 30 40 40 40 40 40];
[10 10 20 20 30 30 30 30 30 40 40 50 60 60 60 60 70];
[10 20 20 20 30 30 40 40 40 50 50 60 70 70 80 80 80];
[20 20 30 30 40 40 50 50 60 60 70 70 80 80 90 100 100];
[20 30 40 40 50 50 60 70 70 80 80 90 100 100 110 120 120];
[30 40 40 50 60 60 70 80 90 90 100 110 110 120 130 140 140];
[40 40 50 60 70 80 90 100 110 120 120 130 140 150 160 170 180];
[40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200];
[50 60 70 80 90 100 120 130 140 150 160 170 180 200 210 220 230];
[50 70 80 90 110 120 130 150 160 170 190 200 210 230 240 260 270];
[60 80 90 100 120 130 150 160 180 190 210 220 240 250 270 280 300];
[70 90 110 120 140 160 170 190 210 230 240 260 280 300 310 330 350];
[80 100 120 130 150 170 190 210 230 250 270 280 300 320 340 360 380];
[80 100 130 150 170 190 210 230 250 270 290 310 330 350 380 400 420];
[90 120 140 160 190 210 230 260 280 310 330 350 380 400 420 440 470];
[100 130 150 180 210 230 250 280 300 330 350 380 400 430 450 480 500];
[110 140 170 200 230 260 290 310 340 370 400 430 460 490 520 540 570];
[120 160 190 220 250 280 310 340 370 410 440 470 500 530 560 590 620];
[140 170 210 240 270 310 340 380 410 440 480 510 550 580 620 650 680];
[150 180 220 250 290 330 360 400 440 470 510 540 580 620 650 690 730];
[150 190 230 270 310 350 380 420 460 490 530 570 610 650 680 720 760];
[160 210 250 290 330 370 410 450 490 530 570 620 660 700 740 780 820];
[180 220 270 310 360 400 440 490 530 580 620 670 710 750 800 840 890];
[180 230 280 320 370 410 460 510 550 600 640 690 740 780 830 880 920];
[200 240 290 340 390 440 490 540 590 640 690 730 780 830 880 930 980];
[200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000];
[220 270 330 380 440 490 550 600 660 710 770 820 880 930 980 1040 1090];
[230 290 350 400 460 520 580 630 690 750 810 860 920 980 1040 1100 1150];
[260 320 390 450 510 580 640 710 770 840 900 960 1030 1090 1160 1220 1290];
[270 330 400 470 540 600 670 730 800 870 930 1000 1070 1130 1200 1260 1330];
[280 350 420 490 560 630 700 770 840 910 980 1050 1120 1190 1260 1330 1400];
[300 380 450 530 600 680 750 830 900 980 1050 1130 1200 1280 1350 1420 1500];
[320 390 480 560 640 720 790 870 950 1030 1110 1190 1270 1350 1430 1510 1590];
[330 420 500 590 670 760 840 920 1010 1090 1170 1260 1340 1430 1510 1600 1680];
[350 430 520 610 700 790 870 960 1050 1130 1220 1310 1400 1480 1570 1660 1740];
[370 460 560 650 740 830 930 1020 1110 1200 1290 1390 1480 1570 1660 1760 1850];
[380 470 570 660 760 850 950 1040 1140 1230 1330 1430 1520 1610 1710 1800 1900];
[390 490 590 690 790 890 990 1090 1190 1290 1390 1490 1590 1690 1780 1880 1980];
[410 520 620 720 830 930 1040 1140 1240 1340 1450 1550 1660 1760 1860 1960 2070];
[430 540 650 760 860 970 1080 1190 1300 1400 1510 1620 1730 1840 1940 2050 2160];
[450 560 670 790 900 1010 1120 1240 1350 1460 1570 1680 1800 1910 2020 2140 2250];
[470 580 700 820 940 1050 1170 1290 1400 1520 1640 1750 1870 1990 2110 2220 2340]],5:50,4:20)

function scribner_volume(small_end_diam,length;decimal_C=false)
    if decimal_C==false
    log_bf=((0.79*small_end_diam^2)-(2*small_end_diam)-4)*(length/16)
    return log_bf
    elseif decimal_C==true
        # ROUND DIB TO NEAREST INTEGER, ROUND LOGFEET DOWN TO NEXT INTEGER
        scale_dib=ceil(Int64,small_end_diam)
        log_feet=floor(Int64,length)
        println(scale_dib,log_feet)
        if scale_dib > 5 && scale_dib < 50 && log_feet > 4 && log_feet < 20
            #STANDARD SIZE LOG, USE LOOKUP TABLE
            log_bf=scribner_decimal_c[scale_dib,log_feet]*1.0
        else
        #OVERSIZE LOG, I.E., DIB>50 AND/OR LOGLENGTH>20, USE FORMULA VERSION
            log_bf=scribner_volume(scale_dib,log_feet,decimal_C=false)
        return log_bf
        end
    end
 end

 #Function to calculate the Doyle scale volume
 #for more info see: http://oak.snr.missouri.edu/forestry_functions/doylebfvolume.php
 function doyle_volume(small_end_diam,length)
   volume=((small_end_diam-4.0)/4.0)^2*length
  return volume
   end

#Function to calculate the International scale volume
#for more info see: http://oak.snr.missouri.edu/forestry_functions/int14bfvolume.php
function international_volume(small_end_diam,length)
  if length == 4
      volume = 0.22*small_end_diam^2-0.71*small_end_diam
    elseif length == 8
      volume = 0.44*small_end_diam^2-1.20*small_end_diam-0.30
    elseif length == 12
      volume = 0.66*small_end_diam^2 - 1.47 * small_end_diam - 0.79
    elseif length == 16
      volume = 0.88 * small_end_diam^2 - 1.52 * small_end_diam - 1.36
    elseif length == 20
      volume = 1.10 * small_end_diam^2 - 1.35 * small_end_diam - 1.90
    elseif length == 24
      volume = 1.10 * small_end_diam^2 - 1.35 * small_end_diam - 1.90 + 0.22 * small_end_diam^2 - 0.71 * small_end_diam
    elseif length == 28
      volume = 1.10 * small_end_diam^2 - 1.35 * small_end_diam - 1.90 + 0.44 * small_end_diam^2 - 1.20 * small_end_diam - 0.30
    elseif length == 32
      volume = 1.10 * small_end_diam^2 - 1.35 * small_end_diam - 1.90 + 0.66 * small_end_diam^2 - 1.47 * small_end_diam - 0.79
    elseif length == 36
      volume = 1.10 * small_end_diam^2 - 1.35 * small_end_diam - 1.90 + 0.88 * small_end_diam^2 - 1.52 * small_end_diam - 1.36
    elseif length == 40
      volume = (1.10 * small_end_diam^2 - 1.35 * small_end_diam - 1.90 )*2
      return volume
    end
  end

####
#I don't know what a NVEL/Volume Equations API should look like?
#Various attempts below...
#introduce abstract super types
abstract type VolumeEquation end
abstract type MerchSpecs end

mutable struct Sawtimber<:MerchSpecs
std_length
trim
min_length
max_length
min_dib
stumpht
end
s=Sawtimber(16.0,0.5,8.0,20.0,6.0,1.0)

mutable struct Fiber<:MerchSpecs
min_length
std_length
min_dib
min_dbh
end
#Fiber(8.0,20.0,2.6,4.6)

mutable struct Pulp<:MerchSpecs
min_live_dib #pulp sawlogs
min_dead_dib #pulp sawlogs
end
#Pulp(6.0,6.0)

####work in progress
#bark thickness separate?
function dib_to_dob() end

#PREDICT DIAMETER INSIDE BARK AT ANY LOGICAL HEIGHT
function get_dib(species,dbh,tht,ht)
end

#PREDICT THE HEIGHT AT WHICH A SPECIFIED DIAMETER INSIDE BARK OCCURS
function get_height_to_dib(species,dbh,tht,dib)
 end

#TO BREAK A SAWTIMBER BOLE INTO LOGS AND POPULATE THE SEGMENT LENGTH
function stem_buck(m::Sawtimber)
    #DETERMINE LENGTH OF BOLE TO DIVIDE INTO LOGS
    bole_length=52.9815
    #bole_length=ht2sawdib-stumpht
    #DETERMINE THE NUMBER OF FULL STANDARD LENGTH LOGS POSSIBLE
    n_logs=floor(Int64,bole_length/(m.std_length+m.trim))
    #POPULATE SEGLEN WITH STANDARD LOG LENGTHS TO START
    seg_length=zeros(n_logs)
    for i in 1:n_logs
        seg_length[i]=m.std_length+m.trim
    end
#ADJUST LOG LENGTHS ACCORDING TO BUCKING RULE AND AMOUNT OF SAWTIMBER EXCESS
    saw_excess=bole_length-n_logs*(m.std_length+m.trim)
    isaw_excess=floor(Int64,saw_excess/2.0)*2
    if n_logs == floor(Int64,n_logs/2.0)*2
        even=true
    else
        even=false
    end
    if isaw_excess >=2.0 && isaw_excess<6.0       #EXTRA SAWEXCESS IS  2' OR 4'
        if even
            if isaw_excess == 2.0
                 seg_length[n_logs-1]=seg_length[n_logs-1]+2.0 # Add 2' to last log
             else
                 seg_length[n_logs-1]=seg_length[n_logs-1]+2.0 # Add 2' to top 2 logs
                 seg_length[n_logs]=seg_length[n_logs]+2.0
             end
        else
            seg_length[n_logs]=seg_length[n_logs] + isaw_excess #Add 2' or 4' to top log
         end
    elseif isaw_excess>=6.0
        if (saw_excess-m.trim) >= m.min_length
            n_logs=n_logs+1 # Create whole new log
            push!(seg_length,isaw_excess+m.trim)
            #seg_length[n_logs]=isaw_excess+m.trim
        else
            n_logs=n_logs+1 # Create new log, split top
            saw_excess=saw_excess+seg_length[n_logs-1]
            seg_length[n_logs]=floor(Int64,saw_excess-(2.0*m.trim)/4.0)*2.0+m.trim
            seg_length[n_logs-1]=floor(saw_excess-seg_length[n_logs]-m.trim/2.0)*2.0+m.trim
        end
    end
    #FAILSAFE TO ENSURE LOGS ARE OF SUFFICIENT LENGTH (B/T MIN SCALE LEN AND MAXLEN)
    sawtop=m.stumpht
    for i in 1:n_logs
        if seg_length[i]<(m.min_length+m.trim)
            seg_length[i]=0
        elseif seg_length[i] > (m.max_length+m.trim)
            seg_length[i]=(m.max_length+m.trim)
        end
        #DETERMINE SAWTOP
        sawtop=seg_length[i]
    return seg_length
    end
end
stem_buck(s)

function stem_buck(m::Fiber)
    fiber_bole=30
#DETERMINE THE NUMBER OF FULL STANDARD LENGTH FIBER LOGS POSSIBLE
n_logs_fiber=floor(Int64,bole_length/m.std_length)

#POPULATE SEGLEN WITH STANDARD LOG LENGTHS TO START
for i in n_logs+1:n_logs+n_logs_fiber
    seg_length[i]=m.std_length
end

#SEE IF ADDITIONAL FIBER ABOVE STANDARD LENGTH LOGS CONSTITUTES A FIBER LOG
fiber_excess=fiber_bole-(n_logs_fiber*m.std_length)
if fiber_excess>m.min_length
    n_logs_fiber=n_logs_fiber+1
    seg_length[n_logs+n_logs_fiber] = fiber_excess
end
end
#TO ASSIGN PRODUCT CLASSES TO EACH LOG SEGMENT IN A SINGLE TREE
function prod_classify()
end

#TO MERCHANDIZE AN INDIVIDUAL TREE WITH PREVIOUSLY BUCKED LOGS AND
#PRODUCTS ASSIGNED TO INDIVIDUAL LOGS INTO PIECES, AND DEVELOP THE NECESSARY PIECE INFORMATION (DIB,LEN,GRS,NET,NPCS)
function merchandize_tree() end


text:
using Dorayaki
using Random
using Plots

const lr = 0.2
const iters = 10000

function main(x, y)
    I, H, O = 1, 10, 1
    l1 = Linear(I, H)
    l2 = Linear(H, O)

    function predict(p)
        q = l1(p)
        q = sigmoid(q)
        q = l2(q)
        return q
    end

    for i in 1:iters
        y_pred = predict(x)
        loss = mean_squared_error(y, y_pred)

        cleargrad!(l1, l2)

        gradient!(loss)

        for l in (l1, l2)
            for p in params(l)
                p.data .-= lr .* p.grad.data
            end
        end

        if i % 1000 == 0
            println(loss.data)
        end
    end
    return predict
end

x = rand(1, 100)
y = sin.(2 * pi .* x) .+ rand(1, 100)

ret = @time main(x, y)
fig = plot(x', y', st=:scatter)
f(x) = ret(x).data[1]
plot!(fig, f)


text:
# Dependancies are called below:
using DataFrames
using CSV
using Statistics

# LEV CREATE
# Ok so here we create the Lev Matrix
# Inputs are two columns that are the two sequences
# Output is the lev Matrix

function lev_create(x,y)

    # ok so we need to get the lengths of both arrays
        lenx = size(x)[1]
        leny = size(y)[1]

    # Initialise the lev array
        lev_array = zeros(Int, leny + 1, lenx + 1)

    #  and now we populate it with the col starting vals
        lev_array[1,:] = 0:(size(lev_array)[2]-1)
        lev_array[:,1] = 0:(size(lev_array)[1]-1)

    # cool so now we do the base lev algo

        for i in 2:size(lev_array)[2]
            for j in 2:size(lev_array)[1]
    #             because we start from 2 not 1, we can skip the first lev condition
                if x[i-1] == y[j-1]
                    sub_cost = 0
                    else sub_cost = 1
                end

                ins = lev_array[j-1,i]+1
                del = lev_array[j,i-1]+1
                sub = lev_array[j-1,i-1] + sub_cost

                lev_array[j,i] = min(ins,del,sub)

            end
        end

        return   lev_array
    end;

# LEV DISTANCE
# This will get the last value of the lev array
# input is a Lev Matrix (or any matrix really) and returns the last value.
# Simple, but well named

function lev_distance(x)

    lev_distance = x[end,end]

        return lev_distance

end

# LEV TRACE
# get a list of the elements that are increasing the edit distance
# This takes a lev matrix, the sequences that created it, and returns an array with
# the edit distance information

# Important to note that this assumes that the master sequence is on the top/x axis
function lev_trace(lev_matrix)

    #     ok so init the bad elements Matrix
        bad_elements = zeros(Int8,1,size(lev_matrix)[2])

    #     ok so start at the bottom right
        i = size(lev_matrix)[2]
        j = size(lev_matrix)[1]
    #     count = 0

         while i > 1 && j > 1

            i_old = i
            j_old = j

            left = lev_matrix[j,i-1]
            diag = lev_matrix[j-1,i-1]
            up   = lev_matrix[j-1,i]

    #         print(left,diag,up)
    #         print(" ")

             if up == min(left,diag,up)
                j = j - 1
    #             print(" up ")
             elseif diag == min(left,diag,up)
                i = i - 1
                j = j - 1
    #             print(" diag ")
            elseif left == min(left,diag,up)
                bad_elements[i]=1
                i = i - 1
    #             print(" left ")

            end
    #             print(i,j)
    #             print(" ")

    #     count = count + 1
        end

    return bad_elements
end

# LEV UPDATE
# This function takes in a lev matrix, slices it at a point and recalculates it from that point
# inputs are: lev matrix, the new sequences, and the positions on the old master sequence were the new sequence changes
# Types are: array, array, array, Int8
function lev_update(matrix,seq,new_master_seq,position)

    #     ok so check how many elements have been removed
        elements_removed = size(matrix)[2] - 1 - size(new_master_seq)[1]

        #     trim this number of elements off the end
        lev_array=matrix[1:end, 1:(size(matrix))[2]-elements_removed]

            for i in (position+1):size(lev_array)[2]
                for j in 2:size(lev_array)[1]

                    if new_master_seq[i-1] == seq[j-1]
                        sub_cost = 0
                        else sub_cost = 1
                    end

                    ins = lev_array[j-1,i]+1
                    del = lev_array[j,i-1]+1
                    sub = lev_array[j-1,i-1] + sub_cost

                    lev_array[j,i] = min(ins,del,sub)
                end
            end

        return lev_array

    end

# BASIC lEV
# This is the basic lev reduction. It takes a Dataframe with the following collumns(order not needed):
# Sequence Identifier
# Sequence Event
# Sequence Order

# Inputs are:
# dataframe, sequence_id, sequence_event, sequence_order
# Types:
# Dataframe, String, String, String

"""
basic_lev(dataframe, sequence_id="sequence_id", sequence_event="sequence_event", sequence_order="sequence_order")

Takes a dataframe of sequences,(type of DataFrames) that has the following three columns:
#sequence_id - the sequence identifier, this might be client, game, or transaction.
#sequence_event - the event or data that is tied to the id, this might be 'new client', 'goal scored' or 'transaction approved'
#sequence_order - the order that the events happened in. this can be any ascending real number.

"""
function basic_lev(dataframe; sequence_id="sequence_id",sequence_event="sequence_event",sequence_order="sequence_order")

    # The lev_master Array is the key object of this funtion. consider it the spine, or memory
    # the Elements are: The Sequence, Lev matrix, Lev Value, Trace
    lev_master = Array[]
    # Elements are: Lev average, and trace average
    lev_history = Array[]

    # ok so first thing is to create a sub dataframe from the dataframe
    gd = groupby(dataframe, Symbol(sequence_id))

    # We need to initalise the master sequence
    master_df = sort!(unique(dataframe[:, [Symbol(sequence_event), Symbol(sequence_order)]]), [Symbol(sequence_order)])
    master_seq=master_df[!,Symbol(sequence_event)]
    # println("Size of master sequence is ",size(master_seq)[1])

    for subdf in gd
    #     Create the first cut of lev matricies
            seq = subdf[!,Symbol(sequence_event)]
            lev_matrix = lev_create(master_seq,seq)
        # Get the lev values
            lev_value = lev_distance(lev_matrix)
        # and do the trace on the lev matrix
            trace = lev_trace(lev_matrix)

            push!(lev_master,[seq,lev_matrix,lev_value,trace])

    end

    # return lev_master

# Cool so this creates the initial lev master object
# Now we want to begin the ole iteration

# Pull out the values so we can do sum statz
lev_value = Int64[]
trace = Array[]
master_seq_hist = Array[]
push!(master_seq_hist,copy(master_seq))

for i in lev_master;
    push!(lev_value,i[3])
    push!(trace,i[4])
end

# Get the history of the run and add it to the lev hist Array
push!(lev_history,Array{Union{Array,String,Float64}}([mean(lev_value)+1,mean(vcat(trace...),dims=1),trace]))
push!(lev_history,Array{Union{Array,String,Float64}}([mean(lev_value)+1,mean(vcat(trace...),dims=1),trace]))
positions = findall(x -> x == maximum(lev_history[end][2]), lev_history[end][2])

# ok so now iterate over this thing
while (lev_history[end][1] < lev_history[end-1][1])
    # println(lev_history[end][1],"   ",lev_history[end-1][1],)

    to_remove=positions[1][2]-1
    # take positions and remove the last one from the master sequence
    # println(" removing event ",master_seq[to_remove]," position ",to_remove," --- ","Avg Lev:",lev_history[end][1],
    #     " Master Sequence Length:",size(master_seq)[1])

    deleteat!(master_seq,to_remove)
    push!(master_seq_hist,copy(master_seq))

    # init these temp arrays
    lev_value = Int64[]
    trace = Array[]
    # update the lev arrays and calculate new lev value
    for i in lev_master;
        i[2] = lev_update(i[2],i[1],master_seq,to_remove)
        i[3] = lev_distance(i[2])
        i[4] = lev_trace(i[2])
        # and create the stats
        push!(lev_value,i[3])
        push!(trace,i[4])
    end

    # Get the history of the run and add it to the lev hist Array
    push!(lev_history,Array{Union{Array,String,Float64}}([mean(lev_value)+1,mean(vcat(trace...),dims=1),trace]))
    positions = findall(x -> x == maximum(lev_history[end][2]), lev_history[end][2])

end;

if lev_history[end][1] > lev_history[end-1][1]
    # println("Going back one and ending things")
    master_seq = master_seq_hist[end-1]
end

    Lev_return = [master_seq,lev_history,lev_master]

    return Lev_return

end


text:
@testset "Jacobian" begin
    graph = Compgraph()
    c = [1.1, 1.2, 1.3, 1.4, 1.5, 1.6]
    add_lincomb!(graph, :B1, c[1], :I, c[2], :A)
    add_lincomb!(graph, :B2, c[3], :I, c[4], :A)
    add_mult!(graph, :N, :B1, :B2)
    add_lincomb!(graph, :D, c[5], :I, c[6], :A)
    add_ldiv!(graph, :R, :D, :N)
    add_output!(graph, :R)

    x = randn(5)
    @test eval_graph(graph, x) ≈
          (c[1] .+ c[2] * x) .* (c[3] .+ c[4] * x) ./ (c[5] .+ c[6] * x)

    dg_dc3 = eval_jac(graph, x, (:B2, 1))
    @test dg_dc3 ≈ (c[1] .+ c[2] * x) ./ (c[5] .+ c[6] * x)

    dg_dc6 = eval_jac(graph, x, (:D, 2))
    @test dg_dc6 ≈
          -x .* (c[1] .+ c[2] * x) .* (c[3] .+ c[4] * x) ./
          (c[5] .+ c[6] * x) .^ 2

    dg_dc2c5 = eval_jac(graph, x, [(:B1, 2), (:D, 1)])
    @test dg_dc2c5 ≈ hcat(
        x .* (c[3] .+ c[4] * x) ./ (c[5] .+ c[6] * x),
        -(c[1] .+ c[2] * x) .* (c[3] .+ c[4] * x) ./ (c[5] .+ c[6] * x) .^ 2,
    )
end


text:
#########
#Created Date: Wednesday April 28th 2021
#Author: Domenico Di Gangi,  <digangidomenico@gmail.com>
#-----
#Last Modified: Thursday June 3rd 2021 5:18:32 pm
#Modified By:  Domenico Di Gangi
#-----
#Description:  Load estimates of latent AR dynamics with model and dgp parameters (variance of latent parameters) equal to that estimated in the application to US congress data
#-----
#-----
########

using DrWatson
@quickactivate "ScoreDrivenExponentialRandomGraphs"
using DataFrames
using PyPlot
pygui(true)
using JLD2
using ScoreDrivenERGM
using ScoreDrivenERGM.DynNets
using Statistics
using Logging
using SharedArrays
using LinearAlgebra
using Logging

@elapsed df = collect_results( datadir("sims", "dgpAR_Fil_SS"))
df["modelTag"] = string.(df["model"])

# begin

# limitSample = 10
limitSample = nothing
winsorProp=0
listParUncMethod = ["WHITE-MLE"]

@elapsed subDf = df[df.T .==73,: ]

@elapsed estResRow = eachrow(subDf)[1]
m = listParUncMethod[1]

for estResRow in eachrow(subDf)

    processFlag = false
    if  (estResRow.T in [73])  &  (estResRow.model.scoreScalingType == "FISH_D") & (estResRow.dgpSettings.type == "AR")
            processFlag = true
    end

    if processFlag
        for m in listParUncMethod

            if winsorProp!=0
                m = "$m$winsorProp"
            end

            timeConf = @elapsed avgCover, constInds, errInds, allConfBandsFiltPar, allConfBandsPar, allmvSDUnParEstCovWhite = average_coverages(estResRow, m; limitSample = limitSample, winsorProp=winsorProp)

            coverDict =  (;m, avgCover, constInds, errInds, allConfBandsFiltPar,  allConfBandsPar, allmvSDUnParEstCovWhite) |>DrWatson.ntuple2dict |> DrWatson.tostringdict

            coverDict = merge(DrWatson.tostringdict(estResRow), coverDict)
            delete!(coverDict, "path")
            delete!(coverDict, "modelTag")

            coverDict["dgp"] = coverDict["dgpSettings"]
            delete!(coverDict, "dgpSettings")
            coverDict["S"] = coverDict["nSample"]
            delete!(coverDict, "nSample")

            # loadPath = estResRow.path

            saveName = replace.( savename(coverDict, "jld2"; allowedtypes = (Real, String, Symbol, NamedTuple, Tuple, ScoreDrivenERGM.DynNets.SdErgm) ), r"[\"]" => "") #loadPath[findlast("\\", loadPath)[end]+1:end]

            timeSave = save( datadir("sims", "dgp&FIl_conf", "cover", saveName), coverDict)

            Logging.@info(" errRatio = $(mean(errInds)) , time = $timeConf")

        end
    end
end

# end

avgCover[:,1,:]
errInds

plot(allConfBandsFiltPar[2,:,1,:,3])
plot(estResRow.allfVecT_filt[2,:,1])
estResRow.allvEstSdResPar[end,:]
names(estResRow)

avgCover[2,1,.!errInds[2,:]]
constInds

@elapsed res = estResRow

N = res.N
T = res.T
nErgmPar = DynNets.number_ergm_par(res.model)
quantilesVals = [ [0.975, 0.025]]
nBands = length(quantilesVals)

if isnothing(limitSample)
    nSample = res.nSample
elseif res.nSample >= limitSample
    nSample = limitSample
else
    error("sample size?")
end

avgCover = SharedArray(zeros(nErgmPar, 2, nSample))
allmvSDUnParEstCovWhite = SharedArray(zeros(3*nErgmPar, 3*nErgmPar,nSample))
allConfBandsPar = SharedArray(zeros(nErgmPar, T, nBands, 2, nSample))
allConfBandsFiltPar = SharedArray(zeros(nErgmPar, T, nBands, 2, nSample))
constInds = SharedArray{Bool}((2, nSample))
errInds = SharedArray{Bool}((2, nSample))

errInds .= true

count = SharedArray(ones(1))

mvSDUnParEstCov = zeros(3,3)
import ScoreDrivenERGM.DynNets:unrestrict_all_par, conf_bands_given_SD_estimates, conf_bands_coverage

Threads.@threads for n=1:nSample

    Logging.@info("Estimating Conf Bands  N = $N , T=$T, $(DynNets.name(res.model)), $(res.dgpSettings) iter n $(count[1]), ")
    count[1] += 1

    vEstSdUnPar = unrestrict_all_par(res.model, res.allvEstSdResPar[:,n])

    ~, allConfBandsFiltPar[:,:,:,:,n], allConfBandsPar[:,:,:,:,n], errFlag, allmvSDUnParEstCovWhite[:,:,n], distribFilteredSD = DynNets.conf_bands_given_SD_estimates(res.model, N, res.allObsT[n], vEstSdUnPar, res.allftot_0[:,n], quantilesVals;  parUncMethod = m, mvSDUnParEstCov=mvSDUnParEstCov, winsorProp=winsorProp )

    coverFiltParUnc = DynNets.conf_bands_coverage(res.allParDgpT[:,:,n],   allConfBandsFiltPar[:,:,:,:,n])

    coverParUnc = DynNets.conf_bands_coverage(res.allParDgpT[:,:,n],   allConfBandsPar[:,:,:,:,n])

    constInds[:, n]  .= any(res.allvEstSdResPar[3:3:end,n].<=0.02, dims=1)
    for indPar in 1:2
        avgCover[indPar, 1, n] =  mean(coverFiltParUnc[indPar,:,1])

        avgCover[indPar, 2, n] =  mean(coverParUnc[indPar,:,1])
    end
        errInds[:, n]  .= errFlag

    # catch
    #     Logging.@warn("Error in estimating confidence bands $T, $N, $n")
    # end

end

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 93c3b8d4-b33e-11eb-335d-2b51a80e19df
using FFTW, Plots, PlutoUI, Random, SampledSignals

# ╔═╡ 0b19c3d5-a429-4fd8-9636-fca0417139f9
include("definitions/SoundPlots.jl")

# ╔═╡ d2818fb6-fb0d-4362-9aed-151f2076a8de
fs = 48000

# ╔═╡ 55962503-6a3f-4caf-99e1-0f3b542ef5e4
md"# White noise
Generating sound using random samples gives white noise.
It has no distinct pitch"

# ╔═╡ 76f78ca5-b6bc-4fb0-9a62-edb03a0b2f17
white = SampleBuf([2*rand()-1 for _ in 1:fs], fs)

# ╔═╡ c0441dfb-58bb-468d-8c41-ca3d514322bf
plotfreqs(white, "White Noise Spectrum")

# ╔═╡ 67790146-df44-40ec-abaf-1b4dfccb8b7b
md"# Periodic Noise
What happens if random noise is repeated?"

# ╔═╡ 30289b86-dd11-45f3-a61b-7257b6539fa9
white1hz = repeat(white, 5)

# ╔═╡ d297644a-27e1-4581-912f-e91f2ee28685
plotfreqs(white1hz, "1 Hz repeated noise")

# ╔═╡ 025f27c1-9dc6-4e7b-ae73-72a298ac93cc
md"With a period of one second the sound is still noise.
You may be able to hear artifacts at regular intervals.
I've listened to several random versions
and it can be easier or harder to pick out the repetition
depending on the particular sound generated.

What about shorter periods?
The repetitions will become faster and more noticable.
At around 30hz it will start to have an audible pitch.
At higher frequencies the pulses smooth out into a buzzy noise with a definite pitch.
At periods of low hundreds of frames the frequency graph shows clearly
the lines of the individual harmonics.
First at the frequency of the repetition,
and with overtones at integer multiples of that frequency."

# ╔═╡ fa4d47c8-6bea-4a57-a9e6-6809c3507b9e
@bind period Slider(100:500:8100)

# ╔═╡ ee695155-4cc7-45b6-8fae-fba482de61b1
period

# ╔═╡ f31f99a8-292e-4525-a6a6-8d0d39a5b088
frequency = fs/period

# ╔═╡ 67007aff-ffd1-4602-8768-80facb50221b
# The sound has a distinct pitch
begin
   reps = fs ÷ period
   base = white[1:period]
   repeated = repeat(base, reps)
end

# ╔═╡ 9b0fa656-52bc-45b8-8721-52ff8d01eb48
plotfreqs(repeated, "Periodic Noise")

# ╔═╡ 9fb202c3-620b-4490-9bfb-25feeda655ae
md"There is no fall off in harmonics,
their amplitude seems randomly distributed arround a flat value.
I was suprised to see how clean the spectrum is
despite using white noise as the base of the sound.
The frequency components of any periodic sound are integer multiples
of the fundamental frequency, which is the inverse of the period."

# ╔═╡ f642114c-5d34-4cf9-843d-35999d2f1587
md"# Odd harmonics with periodic noise
We can remove the even harmonics by making the second half of the waveform
the same as the first half but flipped upside down, call this an odd sound for short. Every even partial is an even repetition, so it splits into two equal parts on the first and second halves of the waveform. To be equal to its own negative on the second half an even partial has to be $0$. Every odd partial reaches the midpoint of the waveform halfway through its own period, and so the second half is the first half offset in phase by half a period. But $\text{sin}(x+\pi) = -\text{sin}(x)$ and so odd harmonics always give odd sounds.
Since a periodic sound is the sum of its even and odd partials
and the odd partials already give an odd sound,
a sound is odd whenever the sum of its even partials are also odd,
but that only happens when the even partials are $0$.

This becomes very useful when imitating natural instruments that only contain odd harmonics, like many wind instruments."

# ╔═╡ 5c8e0d8e-b0a0-4fb6-b43f-41354a4e1ddc
oddnoise = repeat([base; -base], reps÷2)

# ╔═╡ 8c2520cc-e6ad-48d3-bd94-6cb1b8999fed
plotfreqs(oddnoise, "Odd Noise Spectrum")

# ╔═╡ Cell order:
# ╠═93c3b8d4-b33e-11eb-335d-2b51a80e19df
# ╠═d2818fb6-fb0d-4362-9aed-151f2076a8de
# ╠═0b19c3d5-a429-4fd8-9636-fca0417139f9
# ╟─55962503-6a3f-4caf-99e1-0f3b542ef5e4
# ╠═76f78ca5-b6bc-4fb0-9a62-edb03a0b2f17
# ╠═c0441dfb-58bb-468d-8c41-ca3d514322bf
# ╟─67790146-df44-40ec-abaf-1b4dfccb8b7b
# ╠═30289b86-dd11-45f3-a61b-7257b6539fa9
# ╠═d297644a-27e1-4581-912f-e91f2ee28685
# ╟─025f27c1-9dc6-4e7b-ae73-72a298ac93cc
# ╠═fa4d47c8-6bea-4a57-a9e6-6809c3507b9e
# ╠═ee695155-4cc7-45b6-8fae-fba482de61b1
# ╠═f31f99a8-292e-4525-a6a6-8d0d39a5b088
# ╟─67007aff-ffd1-4602-8768-80facb50221b
# ╠═9b0fa656-52bc-45b8-8721-52ff8d01eb48
# ╟─9fb202c3-620b-4490-9bfb-25feeda655ae
# ╟─f642114c-5d34-4cf9-843d-35999d2f1587
# ╠═5c8e0d8e-b0a0-4fb6-b43f-41354a4e1ddc
# ╠═8c2520cc-e6ad-48d3-bd94-6cb1b8999fed

max_stars_repo_name: jnha/music

max_issues_repo_name: jnha/music

max_forks_repo_name: jnha/music


text:
using Documenter, DocumenterLaTeX
using GeoEfficiency

const ROOT = joinpath(@__DIR__, "..")
const PAGES = Any[
    "Home" => "index.md",
    "Manual" => [
		"manual/GeoEfficiency.md",
        #"manual/Error.md",
        #"manual/Input_Console.md",
        "manual/Physics_Model.md",
        "manual/Calculations.md",
        "manual/Output_Interface.md",
        "manual/Input_Batch.md",
     ],

    "Development" => "manual/Development.md",
    "Index" => "list.md",
    "LICENSE" => "LICENSE.md",
]

# GeoEfficiency package docs
doc = makedocs(
    debug 	= true,
    root 	= ROOT,
    build 	= "pdf/build",
    modules = [GeoEfficiency],
    clean 	= "clean" in ARGS,
    format 	= LaTeX(platform = "docker"),
    sitename= "GeoEfficiency.jl",
    authors = "Mohamed E. Krar",
    pages 	= PAGES,
	doctest = "doctest" in ARGS,
    #assets  = ["assets/custom.css"],
);

# hack to only deploy the actual pdf-file
mkpath(joinpath(ROOT, "pdf", "build", "pdfdir"))
let files = readdir(joinpath(ROOT, "pdf", "build"))
    for f in files
        if startswith(f, "GeoEfficiency.jl") && endswith(f, ".pdf")
            mv(joinpath(ROOT, "pdf", "build", f),
               joinpath(ROOT, "pdf", "build", "pdfdir", f))
        end
    end
end

@info "Deploying PDF"
deploydocs(
    repo = "github.com/DrKrar/GeoEfficiency.jl.git",
    root = ROOT,
    target = "pdf/build/pdfdir",
    branch = "gh-pages-pdf",
    forcepush = true,
)


text:
"""
	eitheror : example from Bazinga.jl

	Mathematical program with either-or constraints, from [KMS18].

	Original formulations:

	minimize      (x1 - 8)^2 + (x2 + 3)^2
	subject to    x1 - 2 x2 + 4 ≤ 0    or    x1 - 2 ≤ 0
	              x1^2 - 4 x2 ≤ 0      or    (x1 - 3)^2 + (x2 - 1)^2 - 10 ≤ 0

	Reformulation as a constrained structured problem in the form

	minimize     f(x)
	subject to   c(x) in S

	where
	    f(x) = (x1 - 8)^2 + (x2 + 3)^2
	    c(x) = [ x1 - 2 x2 + 4               ]
       	       [ x1^2 - 4 x2                 ]
   	           [ x1 - 2                      ]
       	       [(x1 - 3)^2 + (x2 - 1)^2 - 10 ]
	    S    = { (a,b) | a ≤ 0  or  b ≤ 0 }^2

    References:
    [KMS18]     Kanzow, Mehlitz, Steck, "Relaxation schemes for mathematical
                programs with switching constraints" (2019).
                ....
                arXiv:1809.02388v1 [math.OC] 7 Sep 2018
"""

using Bazinga, OptiMo
using Random, LinearAlgebra
using DataFrames, CSV
using Printf, Plots

###################################################################################
# problem definition
###################################################################################
mutable struct EITHEROR <: AbstractOptiModel
    meta::OptiModelMeta
end

function EITHEROR()
    name = "eitheror"
    meta = OptiModelMeta(2, 4, x0 = zeros(Float64, 2), name = name)
    return EITHEROR(meta)
end

# necessary methods:
# obj, grad!: cons!, jprod!, jtprod!, proj!, prox!, objprox!
function OptiMo.obj(prob::EITHEROR, x::AbstractVector)
    OptiMo.@lencheck prob.meta.nvar x
    return (x[1] - 8)^2 + (x[2] + 3)^2
end

function OptiMo.grad!(prob::EITHEROR, x::AbstractVector, dfx::AbstractVector)
    OptiMo.@lencheck prob.meta.nvar x
    dfx .= [2 * (x[1] - 8); 2 * (x[2] + 3)]
    return nothing
end

function OptiMo.cons!(prob::EITHEROR, x::AbstractVector, cx::AbstractVector)
    OptiMo.@lencheck prob.meta.nvar x
    OptiMo.@lencheck prob.meta.ncon cx
    cx .= [
        x[1] - 2 * x[2] + 4
        x[1]^2 - 4 * x[2]
        x[1] - 2
        (x[1] - 3)^2 + (x[2] - 1)^2 - 10
    ]
    return nothing
end

function OptiMo.jprod!(
    prob::EITHEROR,
    x::AbstractVector,
    v::AbstractVector,
    Jv::AbstractVector,
)
    OptiMo.@lencheck prob.meta.nvar x v
    OptiMo.@lencheck prob.meta.ncon Jv
    Jv .= [
        v[1] - 2 * v[2]
        2 * x[1] * v[1] - 4 * v[2]
        v[1]
        2 * (x[1] - 3) * v[1] + 2 * (x[2] - 1) * v[2]
    ]
    return nothing
end

function OptiMo.jtprod!(
    prob::EITHEROR,
    x::AbstractVector,
    v::AbstractVector,
    Jtv::AbstractVector,
)
    OptiMo.@lencheck prob.meta.nvar x Jtv
    OptiMo.@lencheck prob.meta.ncon v
    Jtv .= [
        v[1] + 2 * x[1] * v[2] + v[3] + 2 * (x[1] - 3) * v[4]
        -2 * v[1] - 4 * v[2] + 2 * (x[2] - 1) * v[4]
    ]
    return nothing
end

function OptiMo.prox!(prob::EITHEROR, x::AbstractVector, a::Real, z::AbstractVector)
    OptiMo.@lencheck prob.meta.nvar x z
    z .= x
    return nothing
end

function OptiMo.objprox!(prob::EITHEROR, x::AbstractVector, a::Real, z::AbstractVector)
    OptiMo.@lencheck prob.meta.nvar x z
    z .= x
    return 0.0
end

function OptiMo.proj!(prob::EITHEROR, cx::AbstractVector, px::AbstractVector)
    OptiMo.@lencheck prob.meta.ncon cx px
    if cx[1] > 0 && cx[3] > 0
        if cx[1] > cx[3]
            px[1] = cx[1]
            px[3] = 0
        else
            px[1] = 0
            px[3] = cx[3]
        end
    else
        px[1] = cx[1]
        px[3] = cx[3]
    end
    if cx[2] > 0 && cx[4] > 0
        if cx[2] > cx[4]
            px[2] = cx[2]
            px[4] = 0
        else
            px[2] = 0
            px[4] = cx[4]
        end
    else
        px[2] = cx[2]
        px[4] = cx[4]
    end
    #if cx[1] > 0 && cx[3] > 0
    # ...
    return nothing
end

# problem build
problem = EITHEROR()

foldername = "/home/alberto/Documents/Bazinga.jl/demo/data/"
filename = problem.meta.name * "_grid"

# solver build
solver =
    Bazinga.ALPX(max_iter = 50, max_sub_iter = 1000, verbose = false, subsolver = :zerofpr)

R = eltype(problem.meta.x0)
nvar = problem.meta.nvar

xmin = -4.0
xmax = 8.0

data = DataFrame()

xgrid = [(i, j) for i = xmin:0.25:xmax, j = xmin:0.25:xmax];
xgrid = xgrid[:];
ntests = length(xgrid)

x0 = [xgrid[1][1]; xgrid[1][2]]
out = solver(problem, x0 = x0)

for i = 1:ntests
    x0 = [xgrid[i][1]; xgrid[i][2]]

    out = solver(problem, x0 = x0)

    @printf "."
    if mod(i, 50) == 0
        @printf "\n"
    end

    push!(
        data,
        (
            id = i,
            xi_1 = x0[1],
            xi_2 = x0[2],
            xf_1 = out.x[1],
            xf_2 = out.x[2],
            iter = out.iterations,
            time = out.time,
            solved = out.status == :first_order ? 1 : 0,
        ),
    )

end
@printf "\n"

CSV.write(foldername * filename * ".csv", data, header = false)

################################################################################
tolx = 1e-3

global c22 = 0
global c44 = 0
global cun = 0

pyplot()

pts = [(0.0, 0.0)]
rad = sqrt(10.0)
for i = 1:50
    xx = 2.0 * i / 50
    yy = 1.0 - sqrt(10.0 - (xx - 3)^2)
    append!(pts, [(xx, yy)])
end
append!(pts, [(2, 3)])
append!(pts, [(4, 4)])
for i = 1:50
    xx = 4 + (xmax - 4) * i / 50
    yy = (xx^2) / 4
    append!(pts, [(xx, yy)])
end
append!(pts, [(xmin, xmax)])
for i = 0:50
    xx = xmin * (1 - i / 50)
    yy = (xx^2) / 4
    append!(pts, [(xx, yy)])
end
feasset = Shape(pts)
hplt = plot(feasset, color = plot_color(:grey, 0.4), linewidth = 0, legend = false)
xlims!(xmin, xmax)
ylims!(xmin, xmax)

for i = 1:ntests
    xi = [data[i, 2]; data[i, 3]]
    xf = [data[i, 4]; data[i, 5]]

    if norm(xf - [2.0; -2.0], Inf) <= tolx
        global c22 += 1
        scatter!(
            hplt,
            [xi[1]],
            [xi[2]],
            color = :blue,
            marker = :circle,
            markerstrokewidth = 0,
            legend = false,
        )
    elseif norm(xf - [4.0; 4.0], Inf) <= tolx
        global c44 += 1
        scatter!(
            hplt,
            [xi[1]],
            [xi[2]],
            color = :red,
            marker = :diamond,
            markerstrokewidth = 0,
            legend = false,
        )
    else
        global cun += 1
        @printf "(%6.4f,%6.4f) from (%6.4f,%6.4f)\n" xf[1] xf[2] xi[1] xi[2]
    end
end
savefig(foldername * filename * ".pdf")


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to TLE parser.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# Macros tle_str and tlenc_str
# ============================

@testset "Macros tle_str and tlenc_str" begin
    # Read the SCDs TLE from the file.
    tles_file = read_tle("./SCDs.tle")

    # Read the same TLEs from a string.
    tles_str = tle"""
        SCD 1
        1 22490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9991
        2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338066

        SCD 2
        1 25504U 98060A   18165.15074951  .00000201  00000-0  55356-5 0  9994
        2 25504  24.9961  80.1303 0017060 224.4822 286.6438 14.44043397 37312
        """

    # Read the same TLES from a string with wrong checksums.
    # This should not output any exceptions.
    tles_str_nc = tlenc"""
        SCD 1
        1 22490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9990
        2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338060

        SCD 2
        1 25504U 98060A   18165.15074951  .00000201  00000-0  55356-5 0  9990
        2 25504  24.9961  80.1303 0017060 224.4822 286.6438 14.44043397 37310
        """

    # Compare the TLEs.
    @test length(tles_file) == length(tles_str)
    @test length(tles_file) == length(tles_str_nc)

    for i = 1:length(tles_file)
        for sym in fieldnames(TLE)
            @test getfield(tles_file[i], sym) == getfield(tles_str[i], sym)

            # Skip the comparison of checksums for `tles_str_nc`.
            ( (sym == :checksum_l1) || (sym == :checksum_l2) ) && continue
            @test getfield(tles_file[i], sym) == getfield(tles_str_nc[i], sym)
        end
    end
end

# Function: read_tle_from_string
# ==============================

@testset "Function read_tle_from_string" begin
    # Read the SCD 1 TLE from the file.
    tle_scd1_expected = read_tle("./SCDs.tle")[1]

    # Parse the same TLE using the function `read_tle_from_string`.
    tle_scd1_result = read_tle_from_string(
        "1 22490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9991",
        "2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338066")[1]

    # Compare the TLEs.
    for sym in fieldnames(TLE)
        # Skip name and satellite number.
        ( (sym == :name) || (sym == :sat_num) ) && continue

        @test getfield(tle_scd1_result, sym) == getfield(tle_scd1_expected, sym)
    end
end

# Date conversion
# ===============

@testset "Epoch year conversion" begin
    tles = tlenc"""
    1 22490U 93009B   75165.62596833  .00000225  00000-0  11410-4 0  9990
    2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338060

    1 22490U 93009B   76165.62596833  .00000225  00000-0  11410-4 0  9990
    2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338060
    """

    date_1 = JDtoDate(tles[1].epoch)
    date_2 = JDtoDate(tles[2].epoch)

    @test date_1[1] == 2075
    @test date_2[1] == 1976
end

# Test errors
# ===========

@testset "TLE format errors" begin

    @test_throws ErrorException read_tle_from_string("""
    SCD 1
    1 22490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9991
    2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338066
    3 123
    """)

    # Errors in line 1
    # --------------------------------------------------------------------------

    @test_throws ErrorException read_tle_from_string("""
    SCD 1
    1 022490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9991
    2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338066
    """)

    @test_throws ErrorException read_tle_from_string("""
    SCD 1
    1 022490U 93009B   181656.2596833  .00000225  00000-0  11410-4 0  9991
    2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338066
    """)

    # TODO: Errors in line 2.
end

@testset "Checksum errors" begin
    @test_throws ErrorException read_tle("./SCDs-wrong_checksum.tle")
    @test_throws LoadError @eval tle"""
    SCD 1
    1 22490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9991
    2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338066

    SCD 2
    1 25504U 98060A   18165.15074951  .00000201  00000-0  55356-5 0  9994
    2 25504  24.9961  80.1303 0017060 224.4822 286.6438 14.44043397 37313
    """
    @test_throws ErrorException read_tle_from_string(
        "1 22490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9991",
        "2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338063",
        true)
end

# Test conversion to string
# =========================

@testset "Conversion from TLE to string" begin

    # Conversion to string
    # ==========================================================================

    tles = read_tle("./tles_20200122.tle")
    f    = open("./tles_20200122.tle", "r")

    for i = 1:length(tles)
        stri_name = readline(f, keep = true)
        stri_l1   = readline(f, keep = true)
        stri_l2   = readline(f, keep = false)

        # The TLE for AMOS-4 has a `-` before the first time
        # derivative, even though it is 0. Since this is happening only here,
        # we will skip this case.
        tles[i].sat_num == 39237 && continue

        # The conversion of the exponent signal of BSTAR does not have a
        # defined pattern if BSTAR = 0.
        bstar_exp_le = false
        if tles[i].bstar == 0
            bstar_exp_le = stri_l1[60] == '-'
        end

        stri = stri_name * stri_l1 * stri_l2
        strf = tle_to_str(tles[i], bstar_exp_le = bstar_exp_le)

        @test strf == stri
    end

    close(f)

    # Printing to IO
    # ==========================================================================

    str_tle_cbers4 = """
    CBERS 4
    1 40336U 14079A   18166.15595376 -.00000014  00000-0  10174-4 0  9993
    2 40336  98.4141 237.7928 0001694  75.7582 284.3804 14.35485112184485"""

    tle_cbers4 = read_tle_from_string(str_tle_cbers4)[1]

    tle_wcs = TLE(tle_cbers4; checksum_l1 = 0, checksum_l2 = 1)

    # Print without modifying the checksum.
    io = IOBuffer()
    print_tle(io, tle_wcs, recompute_checksum = false)
    result = String(take!(io))

    expected = """
    CBERS 4
    1 40336U 14079A   18166.15595376 -.00000014  00000-0  10174-4 0  9990
    2 40336  98.4141 237.7928 0001694  75.7582 284.3804 14.35485112184481"""

    @test result == expected

    # Print recomputing the checksum.
    io = IOBuffer()
    print_tle(io, tle_wcs)
    result = String(take!(io))
    expected = str_tle_cbers4

    @test result == expected
end


text:
using Discreet
import StatsBase: FrequencyWeights, ProbabilityWeights
using Compat, Compat.Test
import Compat.sum

freqs_uniform = FrequencyWeights([1, 1, 1, 1, 1, 1])
probs_uniform = ProbabilityWeights([.5, .5])
sample = [:a, :b, :c, :d, :e, :f]

@testset "Entropy routines" begin
    ## Naive estimator
    @test log(6) ≈ entropy(freqs_uniform)
    @test 0 == entropy(FrequencyWeights([6]))
    @test log(6) ≈ estimate_entropy(sample)

    ## CS and shrinkage estimates, from library(entropy) in R
    @test 3.840549310406 ≈ entropy(FrequencyWeights([1,1,1,1,1,1]); method=:CS)
    @test 2.201137101279 ≈ entropy(FrequencyWeights([4,2,3,2,4,2,1,1]); method=:CS)
    @test 1.968382408728 ≈ entropy(FrequencyWeights([4,2,3,2,4,2,1,1]))
    @test 2.379602895309 ≈ entropy(FrequencyWeights([4,2,3,0,2,4,0,0,2,1,1]);
                                   method=:Shrink)

    ## Any other method
    @test_throws ArgumentError entropy(freqs_uniform; method=:NotImplemented)
    @test_throws ArgumentError estimate_entropy(sample; method=:NotImplemented)

    ## Joint entropy
    @test log(6) ≈ estimate_joint_entropy(sample, sample)
    @test 3.840549310406 ≈ estimate_joint_entropy(sample, sample; method=:CS)
end

@testset "Mutual information routines" begin
    labels_a = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]
    labels_b = [1, 1, 1, 1, 2, 1, 2, 2, 2, 2, 3, 1, 3, 3, 3, 2, 2]

    mi = mutual_information(labels_a, labels_b)
    @test 0.41022 ≈ mi atol=1e-5

    e_a = estimate_entropy(labels_a)
    e_b = estimate_entropy(labels_b)
    mi_normalized = mutual_information(labels_a, labels_b; normalize=true)
    @test mi / min(e_a, e_b) ≈ mi_normalized

    data = hcat(sample, sample, sample)
    @test log(6) * ones(3, 3) ≈ mutual_information(data)
end

@testset "Contingency tables" begin
    d = rand(10, 10)
    table = d / sum(d)
    mi_1 = mutual_information_contingency(table)
    e_a, e_b =  entropy(sum(table, dims=1)), entropy(sum(table, dims=2))
    mi_2 = e_a + e_b - entropy(table[:])
    @test mi_1 ≈ mi_2

    ## Normalized MI
    mi_norm = mutual_information_contingency(table; normalize=true)
    @test mi_norm ≈ mi_1 / min(e_a, e_b)

    ## Table with no association
    @test 0 ≈ mutual_information_contingency(ones(20, 5) / 100) atol=1e-10
end


text:
"""
    ErrorMap

A nonlinear mapping between the space of unit quaternions and three-dimensional
rotation errors.

These mappings are extremely useful for converting from globally nonsingular 3D rotation
representations such as `UnitQuaternion` or `RotMatrix3` to a three-parameter error that can
be efficiently used in gradient-based optimization methods that optimize deviations about
a current iterate using first or second-order information.

# Usage
    errmap(v::AbstractVector)  # "forward" map from a 3D error to a `UnitQuaternion`
    inv(errmap)(R::Rotation)   # "inverse" map from a rotation (via `UnitQuaternion`) to a 3D error

where `errmap <: ErrorMap`

# Implemented Maps
- `CayleyMap`:  Uses `RodriguesParam` as the error representation (default).
    Goes singular at 180° and does not have a sign ambiguity.
- `ExponentialMap`: Uses the canonical exponential map from Lie Group theory. Computationally
    expensive to compute. Exhibits kinematic singularities.
- `MRPMap`: Uses a scaled `MRP` as the error representation. Singular at 360° but has
    a sign ambiguity (with the "shadow" MRP set).
- `QuatVecMap`: Uses the vector part of the quaternions. Cheapest map to compute, but
    goes singular at 180° and suffers from sign ambiguity.
- `IdentityMap`: Maps values through directly. Only works with three-parameter rotation
    representations with the following methods: `R(::SVector{3})` and `SVector(::R)::SVector{3}`
"""
abstract type ErrorMap end
struct CayleyMap <: ErrorMap end
struct ExponentialMap <: ErrorMap end
struct MRPMap <: ErrorMap end
struct QuatVecMap <: ErrorMap end
struct IdentityMap <: ErrorMap end

"""
    InvErrorMap

The nonlinear mapping from unit quaternions to a three-dimensional error state. Obtained by
inverting an `ErrorMap`, i.e.

    InvCayleyMap() = inv(CayleyMap())

# Usage
    imap(R::Rotation)             # "inverse" map from a rotation to a 3D error
    inv(imap)(v::AbstractVector)  # "forward" map from a 3D error to a `UnitQuaternion`

where `imap <: InvErrorMap`.

See `ErrorMap` for documentation on the implemented maps.
"""
abstract type InvErrorMap end
struct InvCayleyMap <: InvErrorMap end
struct InvExponentialMap <: InvErrorMap end
struct InvMRPMap <: InvErrorMap end
struct InvQuatVecMap <: InvErrorMap end

inv(::CayleyMap) = InvCayleyMap()
inv(::ExponentialMap) = InvExponentialMap()
inv(::QuatVecMap) = InvQuatVecMap()
inv(::MRPMap) = InvMRPMap()
inv(::IdentityMap) = IdentityMap()

inv(::InvCayleyMap) = CayleyMap()
inv(::InvExponentialMap) = ExponentialMap()
inv(::InvQuatVecMap) = QuatVecMap()
inv(::InvMRPMap) = MRPMap()

# Scalings
@inline scaling(::Type{ExponentialMap}) = 0.5
@inline scaling(::Type{QuatVecMap}) = 1.0
@inline scaling(::Type{CayleyMap}) = 1.0
@inline scaling(::Type{MRPMap}) = 2.0
scaling(m::M) where M <: ErrorMap = scaling(M)

# Quaternion Maps
(::ExponentialMap)(ϕ::AbstractVector) = expm(ϕ/scaling(ExponentialMap))

function (::QuatVecMap)(v::AbstractVector)
    check_length(v, 3)
    μ = 1/scaling(QuatVecMap)
    UnitQuaternion(sqrt(1-μ^2*v'v), μ*v[1], μ*v[2], μ*v[3])
end

function (::CayleyMap)(g::AbstractVector)
    check_length(g, 3)
    g /= scaling(CayleyMap)
    M = 1/sqrt(1+g'g)
    UnitQuaternion(M, M*g[1], M*g[2], M*g[3])
end

function (::MRPMap)(p::AbstractVector)
    check_length(p, 3)
    p /= scaling(MRPMap)
    n2 = p'p
    M = 2/(1+n2)
    UnitQuaternion((1-n2)/(1+n2), M*p[1], M*p[2], M*p[3])
end

function (::IdentityMap)(q::AbstractVector)
    check_length(q, 3)
    UnitQuaternion(q[1], q[2], q[3], q[4])
end

# Quaternion Map Jacobians
"""
    jacobian(::ErrorMap, ϕ)

Jacobian of the quaternion map that takes a three-dimensional vector `ϕ` and returns a
    unit quaternion.
Returns a 4x3 Static Matrix

For all the maps (except the `IdentityMap`)
`jacobian(::ErrorMap, zeros(3)) = [0; I] = Hmat()'`
"""
function jacobian(::ExponentialMap, ϕ, eps=1e-5)
    μ = 1/scaling(ExponentialMap)
    θ = norm(ϕ)
    cθ = cos(μ*θ/2)
    sincθ = sinc(μ*θ/2π)
    if θ < eps
        0.5*μ*[-0.5*μ*sincθ*ϕ'; sincθ*I + (cθ - sincθ)*ϕ*ϕ']
    else
        0.5*μ*[-0.5*μ*sincθ*ϕ'; sincθ*I + (cθ - sincθ)*ϕ*ϕ'/(ϕ'ϕ)]
    end
end

function jacobian(::QuatVecMap, v)
    μ = 1/scaling(QuatVecMap)
    μ2 = μ*μ
    M = -μ2/sqrt(1-μ2*v'v)
    @SMatrix [v[1]*M v[2]*M v[3]*M;
              μ 0 0;
              0 μ 0;
              0 0 μ]
end

function jacobian(::CayleyMap, g)
    μ = 1/scaling(CayleyMap)
    μ2 = μ*μ
    n = 1+μ2*g'g
    ni = 1/n
    μ*[-μ*g'; -μ2*g*g' + I*n]*ni*sqrt(ni)
end

function jacobian(::MRPMap, p)
    μ = 1/scaling(MRPMap)
    μ2 = μ*μ
    n = 1+μ2*p'p
    2*[-2*μ2*p'; I*μ*n - 2*μ*μ2*p*p']/n^2
end

jacobian(::IdentityMap, q) = I

############################################################################################
#                             INVERSE RETRACTION MAPS
############################################################################################
# (emap::ErrorMap)(R::Rotation) =
#     throw(ArgumentError("Must use the inverse map to convert a rotation to an error state, e.g. inv(emap)(q)"))
# (emap::InvErrorMap)(R::Rotation) = emap(UnitQuaternion(R))  # automatically call the inverse map for rotations
# NOTE: Julia v1.0 doesn't allow these methods on abstract types. Leave out for now.

for imap in (InvExponentialMap,InvCayleyMap,InvMRPMap,InvQuatVecMap)
    @eval begin
        @inline (imap::$imap)(R::Rotation) = imap(UnitQuaternion(R))
    end
end

(::InvExponentialMap)(q::UnitQuaternion) = scaling(ExponentialMap)*logm(q)
(::InvCayleyMap)(q::UnitQuaternion) = scaling(CayleyMap) * vector(q)/q.w
(::InvMRPMap)(q::UnitQuaternion) = scaling(MRPMap)*vector(q)/(1+q.w)
(::InvQuatVecMap)(q::UnitQuaternion) = scaling(QuatVecMap)*vector(q) * sign(q.w)

# ~~~~~~~~~~~~~~~ Inverse map Jacobians ~~~~~~~~~~~~~~~ #
"""
    jacobian(::InvErrorMap, q::UnitQuaternion)

Jacobian of the inverse quaternion map, returning a 3x4 matrix.
For all maps: `jacobian(::InvErrorMap, UnitQuaternion(I)) = [0 I] = Hmat()'`
"""
function jacobian(::InvExponentialMap, q::UnitQuaternion, eps=1e-5)
    μ = scaling(ExponentialMap)
    s = scalar(q)
    v = vector(q)
    θ2 = v'v
    θ = sqrt(θ2)
    datan = 1/(θ2 + s^2)
    ds = -datan*v

    if θ < eps
        return 2*μ*[ds (v*v' + I)/s]
    else
        atanθ = atan(θ,s)
        dv = ((s*datan - atanθ/θ)v*v'/θ + atanθ*I )/θ
        d0 = ((s*datan - atanθ/θ)v*v'/θ^2 + atanθ/θ*I )
        d1 = (s*datan - atanθ/θ)
        d2 = v*v'/θ2
        d3 = atanθ/θ * I
        return 2*μ*[ds dv]
    end
end

function jacobian(::InvQuatVecMap, q::UnitQuaternion)
    μ = scaling(QuatVecMap)
    return sign(q.w) * SA[
                0. μ 0 0;
                0. 0 μ 0;
                0. 0 0 μ]
end

function jacobian(::InvCayleyMap, q::UnitQuaternion)
    μ = scaling(CayleyMap)
    si = 1/q.w
    return μ*@SMatrix [-si^2*q.x si 0 0;
                       -si^2*q.y 0 si 0;
                       -si^2*q.z 0 0 si]
end

function jacobian(::InvMRPMap, q::UnitQuaternion)
    μ = scaling(MRPMap)
    si = 1/(1+q.w)
    return μ*@SMatrix [-si^2*q.x si 0 0;
                       -si^2*q.y 0 si 0;
                       -si^2*q.z 0 0 si]
end

jacobian(::IdentityMap, q::UnitQuaternion) = I

# ~~~~~~~~~~~~~~~ Inverse map Jacobian derivative ~~~~~~~~~~~~~~~ #
"""
    ∇jacobian(::InvErrorMap, q::UnitQuaternion, b::SVector{3})

Jacobian of G(q)'b, where G(q) = jacobian(::InvErrorMap, q),
    b is a 3-element vector
"""
function ∇jacobian(::InvExponentialMap, q::UnitQuaternion, b::SVector{3}, eps=1e-5)
    μ = scaling(ExponentialMap)
    s = scalar(q)
    v = vector(q)
    θ2 = v'v
    θ = sqrt(θ2)
    datan = 1/(θ2 + s^2)
    ds = -datan*v

    if θ < eps
        # return 2*μ*[b'ds; (v*v'b + b)/s]
        return 2*μ*[b'*(datan^2*2s*v) -b'datan*I;
                    -(v*v'b +b)/s^2 (I*(v'b) + v*b')/s]
    else
        dsds = 2b's*datan^2*v
        dsdv = b'*(-datan*I + 2datan^2*v*v')

        atanθ = atan(θ,s)
        d1 = (s*datan - atanθ/θ)
        d2 = v*v'b/θ2
        d3 = atanθ/θ*b
        d1ds = (datan - 2s^2*datan^2 + datan)
        dvds = d1ds*d2 - datan*b

        d1dv =  (-2s*datan^2*v' - s*datan*v'/θ^2 + atanθ/θ^3*v')
        d2dv = (I*(v'b) + v*b')/θ2 - 2(v*v'b)/θ^4 * v'
        d3dv = b*(s*datan*v'/θ^2 - atanθ/θ^3*v')
        dvdv = d2*d1dv + d1*d2dv + d3dv

        # return 2*μ*[ds'b; dv'b]
        return 2*μ*@SMatrix [
            dsds    dsdv[1] dsdv[2] dsdv[3];
            dvds[1] dvdv[1] dvdv[4] dvdv[7];
            dvds[2] dvdv[2] dvdv[5] dvdv[8];
            dvds[3] dvdv[3] dvdv[6] dvdv[9];
        ]
    end
end

function ∇jacobian(::InvCayleyMap, q::UnitQuaternion, b::SVector{3})
    μ = scaling(CayleyMap)
    si = 1/q.w
    v = vector(q)
    μ*@SMatrix [
        2*si^3*(v'b) -si^2*b[1] -si^2*b[2] -si^2*b[3];
       -si^2*b[1] 0 0 0;
       -si^2*b[2] 0 0 0;
       -si^2*b[3] 0 0 0;
    ]
end

function ∇jacobian(::InvMRPMap, q::UnitQuaternion, b::SVector{3})
    μ = scaling(MRPMap)
    si = 1/(1+q.w)
    v = vector(q)
    μ * @SMatrix [
        2*si^3*(v'b) -si^2*b[1] -si^2*b[2] -si^2*b[3];
       -si^2*b[1] 0 0 0;
       -si^2*b[2] 0 0 0;
       -si^2*b[3] 0 0 0;
    ]
end

function ∇jacobian(::InvQuatVecMap, q::UnitQuaternion, b::SVector{3})
    μ = scaling(QuatVecMap)
    @SMatrix zeros(4,4)
end


text:
# related to
# https://github.com/JuliaDiffEq/DiffEqParamEstim.jl/blob/72114707e95e6a19dd264c8fdbb476e9fda6ee31/src/multiple_shooting_objective.jl
using KrylovKit, DiffEqBase

abstract type AbstractShootingProblem end

####################################################################################################
function flow(x, tspan, pb::ODEProblem; alg = Euler(), kwargs...)
	_prob = DiffEqBase.remake(pb; u0 = x, tspan = tspan)
	# the use of concrete_solve makes it compatible with Zygote
	sol = DiffEqBase.concrete_solve(_prob, alg, save_everystep = false; kwargs...)[end]
	return sol
end

# this function takes into accound a parameter passed to the vector field
# Putting the options `save_start=false` seems to bug with Sundials
function flow(x, p, tspan, pb::ODEProblem; alg = Euler(), kwargs...)
	_prob = DiffEqBase.remake(pb; u0 = x, tspan = tspan, p = p)
	# the use of concrete_solve makes it compatible with Zygote
	sol = DiffEqBase.concrete_solve(_prob, alg, save_everystep = false; kwargs...)[end]
	return sol
end
################
# this flow is used for computing the derivative of the flow, so pb encode the variational equation
function dflow(dx, x::AbstractVector, tspan, pb::ODEProblem; alg = Euler(), kwargs...)
	n = length(x)
	_prob = DiffEqBase.remake(pb; u0 = vcat(x, dx), tspan = tspan)
	# the use of concrete_solve makes it compatible with Zygote
	sol = DiffEqBase.concrete_solve(_prob, alg, save_everystep = false; kwargs...)[end]
	return sol[1:n], sol[n+1:end]
end

function dflow(dx, x::AbstractVector, p, tspan, pb::ODEProblem; alg = Euler(), kwargs...)
	n = length(x)
	_prob = DiffEqBase.remake(pb; u0 = vcat(x, dx), tspan = tspan, p = p)
	# the use of concrete_solve makes it compatible with Zygote
	sol = DiffEqBase.concrete_solve(_prob, alg, save_everystep = false; kwargs...)[end]
	return sol[1:n], sol[n+1:end]
end
################
function dflow_fd(x, dx, tspan, pb::ODEProblem; alg = Euler(), δ = 1e-9, kwargs...)
	sol1 = flow(x .+ δ .* dx, tspan, pb; alg = alg, kwargs...)
	sol2 = flow(x 			, tspan, pb; alg = alg, kwargs...)
	return sol2, (sol1 .- sol2) ./ δ
end

# this function takes into accound a parameter passed to the vector field
function dflow_fd(x, dx, p, tspan, pb::ODEProblem; alg = Euler(), δ = 1e-9, kwargs...)
	sol1 = flow(x .+ δ .* dx, p, tspan, pb; alg = alg, kwargs...)
	sol2 = flow(x 			, p, tspan, pb; alg = alg, kwargs...)
	return sol2, (sol1 .- sol2) ./ δ
end
################
# this gives access to the full solution, convenient for Poincaré shooting
function flowFull(x, tspan, pb::ODEProblem; alg = Euler(), kwargs...)
	_prob = DiffEqBase.remake(pb; u0 = x, tspan = tspan)
	sol = DiffEqBase.solve(_prob, alg; kwargs...)
end

# this function takes into accound a parameter passed to the vector field and returns the full solution from the ODE solver. This is useful in Poincare Shooting to extract the period.
function flowFull(x, p, tspan, pb::ODEProblem; alg = Euler(), kwargs...)
	_prob = DiffEqBase.remake(pb; u0 = x, tspan = tspan, p = p)
	sol = DiffEqBase.solve(_prob, alg; kwargs...)
end

####################################################################################################
# Structures related to computing ODE/PDE Flows
"""
	fl = Flow(F, flow, dflow)
This composite type encapsulates:
- the vector field `x -> F(x)` associated to a Cauchy problem,
- the flow (or semigroup) associated to the Cauchy problem `(x, t) -> flow(x, t)`. Only the last time point must be returned.
- the flow (or semigroup) associated to the Cauchy problem `(x, t) -> flow(x, t)`. The whole solution on the time interval (0,t) must be returned. This is not strictly necessary to provide this.
- the differential `dflow` of the flow w.r.t. `x`, `(x, dx, t) -> dflow(x, dx, t)`. One important thing is that we require `dflow(x, dx, t)` to return 2 vectors: the first is `flow(x, t)` and the second is the value of the derivative of the flow as the second vector.

There are some simple constructors for which you only have to pass a `prob::ODEProblem` from `DifferentialEquations.jl` and an ODE time stepper like `Tsit5()`. Hence, you can do for example

	fl = Flow(F, prob, Tsit5(); kwargs...)

If your vector field depends on parameters `p`, you can define a `Flow` using

	fl = Flow(F, p, prob, Tsit5(); kwargs...)

Finally, you can pass two `ODEProblem` where the second one is used to compute the variational equation:

	fl = Flow(F, p, prob1::ODEProblem, alg1, prob2::ODEProblem, alg2; kwargs...)

"""
struct Flow{TF, Tf, Tff, Td}
	F::TF				# vector field F(x)
	flow::Tf			# flow(x,t)
	flowFull::Tff		# flow(x,t) returns full solution
	dflow::Td			# dflow(x,dx,t) returns (flow(x,t), dflow(x,t)⋅dx)

	function Flow(F::TF, fl::Tf, flf::Tff, df::Td) where {TF, Tf, Tff, Td}
		new{TF, Tf, Tff, Td}(F, fl, flf, df)
	end

	function Flow(F::TF, fl::Tf, df::Td) where {TF, Tf, Td}
		new{TF, Tf, Nothing, Td}(F, fl, nothing, df)
	end
end

(fl::Flow)(x, tspan)     = fl.flow(x, tspan)
(fl::Flow)(x, dx, tspan) = fl.dflow(x, dx, tspan)

"""
Creates a Flow variable based on a `prob::ODEProblem` and ODE solver `alg`. The vector field `F` has to be passed, this will be resolved in the future as it can be recovered from `prob`. Also, the derivative of the flow is estimated with finite differences.
"""
function Flow(F, prob::ODEProblem, alg; kwargs...)
	return Flow(F,
		(x, t) -> 				flow(x, (zero(t), t), prob; alg = alg, kwargs...),
		(x, t) ->			flowFull(x, (zero(t), t), prob; alg = alg, kwargs...),
		(x, dx, t) -> dflow_fd((x, dx), (zero(t), t), prob; alg = alg, kwargs...)
		)
end

# this constructor takes into accound a parameter passed to the vector field
function Flow(F, p, prob::ODEProblem, alg; kwargs...)
	return Flow(F,
		(x, t) ->			 	flow(x, p, (zero(t), t), prob; alg = alg, kwargs...),
		(x, t) -> 		 	flowFull(x, p, (zero(t), t), prob; alg = alg, kwargs...),
		(x, dx, t) -> dflow_fd((x, dx), p, (zero(t), t), prob; alg = alg, kwargs...)
		)
end

function Flow(F, p, prob1::ODEProblem, alg1, prob2::ODEProblem, alg2; kwargs...)
	return Flow(F,
		(x, t) -> 			flow(x, p, (zero(t), t), prob1, alg = alg1; kwargs...),
		(x, t) -> 		flowFull(x, p, (zero(t), t), prob1, alg = alg1; kwargs...),
		(x, dx, t) ->  dflow(dx, x, p, (zero(t), t), prob2, alg = alg2; kwargs...)
		)
end

"""
Creates a Flow variable based on a `prob1::ODEProblem` and ODE solver `alg`. The derivative of the flow is imlemented in using `prob2::ODEProblem`. The vector field `F` has to be passed, this will be resolved in the future as it can be recovered from `prob`.
"""
function Flow(F, prob::ODEProblem, alg; kwargs...)
	return Flow(F,
		(x, t) ->			  flow(x, (zero(t), t), prob; alg = alg, kwargs...),
		(x, t) ->		  flowFull(x, (zero(t), t), prob; alg = alg, kwargs...),
		(x, dx, t) -> dflow_fd(x, dx, (zero(t), t), prob; alg = alg, kwargs...)
		)
end

# this constructor takes into accound a parameter passed to the vector field
function Flow(F, p, prob::ODEProblem, alg; kwargs...)
	return Flow(F,
		(x, t) ->			  flow(x, p, (zero(t), t), prob; alg = alg, kwargs...),
		(x, t) -> 		  flowFull(x, p, (zero(t), t), prob; alg = alg, kwargs...),
		(x, dx, t) -> dflow_fd(x, dx, p, (zero(t), t), prob; alg = alg, kwargs...)
		)
end
####################################################################################################
# Standard Shooting functional
"""
	pb = ShootingProblem(flow::Flow, ds, section)

This composite type implements the Standard Simple / Multiple Standard Shooting  method to locate periodic orbits. The arguments are as follows
- `flow::Flow`: implements the flow of the Cauchy problem though the structure `Flow`.
- `ds`: vector of time differences for each shooting. Its length is written `M`. If `M==1`, then the simple shooting is implemented and the multiple one otherwise.
- `section`: implements a phase condition. The evaluation `section(x)` must return a scalar number where `x` is a guess for the periodic orbit. Note that the period `T` of the guess `x` is always included either as the last component of `T = x[end]` or as `T = x.p`. The type of `x` depends on what is passed to the newton solver.

You can then call `pb(orbitguess)` to apply the functional to a guess. Note that `orbitguess` must be of size M * N + 1 where N is the number of unknowns of the state space and `orbitguess[M * N + 1]` is an estimate of the period `T` of the limit cycle. This form of guess is convenient for the use of the linear solvers in `IterativeSolvers.jl` (for example) which accepts only `AbstractVector`s. Another accepted guess is of the form `BorderedArray(guess, T)` where `guess[i]` is the state of the orbit at the `i`th time slice. This last form allows for non-vector state space which can be convenient for 2d problems for example.

A functional, hereby called `G`, encodes the shooting problem. For example, the following methods are available:
- `pb(orbitguess)` evaluates the functional G on `orbitguess`
- `pb(orbitguess, du)` evaluates the jacobian `dG(orbitguess).du` functional at `orbitguess` on `du`

## Simplified constructors
- A simpler way to build a functional is to use
	pb = ShootingProblem(F, p, prob::ODEProblem, alg, M::Int, section; kwargs...)
or

	pb = ShootingProblem(F, p, prob::ODEProblem, alg, ds, section; kwargs...)
where `F` is the vector field, `p` is a parameter (to be passed to the vector and the flow), `prob` is an `ODEProblem` which is used to create a flow using the ODE solver `alg` (for example `Tsit5()`). Finally, the arguments `kwargs` are passed to the ODE solver defining the flow. Look at `DifferentialEquations.jl` for more information. Note that, in this case, the derivative of the flow is computed internally using Finite Differences.
- The other way is an elaboration of the previous one
	pb = ShootingProblem(F, p, prob1::ODEProblem, alg1, prob2::ODEProblem, alg2, M::Int, section; kwargs...)
or

	pb = ShootingProblem(F, p, prob1::ODEProblem, alg1, prob2::ODEProblem, alg2, ds, section; kwargs...)
where we supply now two `ODEProblem`. The first one `prob1`, is used to define the flow associated to `F` while the second one is a problem associated to the derivative of the flow. Hence, `prob2` should implement the following vector field ``\\tilde F(x,\\phi) = (F(x),dF(x)\\cdot\\phi)``.
"""
@with_kw struct ShootingProblem{Tf <: Flow, Ts, Tsection} <: AbstractShootingProblem
	flow::Tf							# should be a Flow{TF, Tf, Td}
	ds::Ts = diff(LinRange(0, 1, 5))	# number of Poincaré sections
	section::Tsection					# sections for phase condition
end

ShootingProblem(F, prob::ODEProblem, alg, ds, section; kwargs...) =  ShootingProblem(Flow(F, prob, alg; kwargs...), ds, section)

ShootingProblem(F, prob::ODEProblem, alg, M::Int, section; kwargs...) = ShootingProblem(Flow(F, prob, alg; kwargs...), diff(LinRange(0, 1, M + 1)), section)

# this constructor takes into accound a parameter passed to the vector field
ShootingProblem(F, p, prob::ODEProblem, alg, ds, section; kwargs...) = ShootingProblem(Flow(F, p, prob, alg; kwargs...), ds, section)

ShootingProblem(F, p, prob::ODEProblem, alg, M::Int, section; kwargs...) = ShootingProblem(Flow(F, p, prob, alg; kwargs...), diff(LinRange(0, 1, M + 1)), section)

# idem but with an ODEproblem to define the derivative of the flow
ShootingProblem(F, p, prob1::ODEProblem, alg1, prob2::ODEProblem, alg2, ds, section; kwargs...) = ShootingProblem(Flow(F, p, prob1, alg1, prob2, alg2; kwargs...), ds, section)

ShootingProblem(F, p, prob1::ODEProblem, alg1, prob2::ODEProblem, alg2, M::Int, section; kwargs...) = ShootingProblem(Flow(F, p, prob1, alg1, prob2, alg2; kwargs...), diff(LinRange(0, 1, M + 1)), section)

# this function extracts the last component of the periodic orbit
extractPeriodShooting(x::AbstractVector) = x[end]
extractPeriodShooting(x::BorderedArray)  = x.p

# Standard shooting functional
function (sh::ShootingProblem)(x::AbstractVector)
	# period of the cycle
	# Sundials does not like @views :(
	T = extractPeriodShooting(x)
	M = length(sh.ds)
	N = div(length(x) - 1, M)

	# extract the orbit guess and reshape it into a matrix as it's more convenient to handle it
	xv = @view x[1:end-1]
	xc = reshape(xv, N, M)

	# variable to hold the computed result
	out = similar(x)
	outv = @view out[1:end-1]
	outc = reshape(outv, N, M)

	for ii in 1:M-1
		# yi = Flow(xi, dsi⋅T) - x[i+1]
		outc[:, ii] .= sh.flow(xc[:, ii], sh.ds[ii] * T) .- xc[:, ii+1]
	end
	# ym = Flow(xm, dsm⋅T) - x1
	outc[:, M] .= sh.flow(xc[:, M], sh.ds[M] * T) .- xc[:, 1]

	# add constraint
	out[end] = sh.section(x)

	return out
end

# jacobian of the shooting functional
function (sh::ShootingProblem)(x::AbstractVector, dx::AbstractVector; δ = 1e-8)
	# period of the cycle
	# Sundials does not like @views :(
	dT = extractPeriodShooting(dx)
	T  = extractPeriodShooting(x)
	M = length(sh.ds)
	N = div(length(x) - 1, M)

	xv = @view x[1:end-1]
	xc = reshape(xv, N, M)

	dxv = @view dx[1:end-1]
	dxc = reshape(dxv, N, M)

	# variable to hold the computed result
	out = similar(x)
	outv = @view out[1:end-1]
	outc = reshape(outv, N, M)

	for ii = 1:M-1
		# call jacobian of the flow
		tmp = sh.flow(xc[:, ii], dxc[:, ii], sh.ds[ii] * T)
		outc[:, ii] .= tmp[2] .+ sh.flow.F(tmp[1]) .* sh.ds[ii] * dT .- dxc[:, ii+1]
	end

	# ym = Flow(xm, dsm⋅T) - x1
	tmp = sh.flow(xc[:, M], dxc[:, M], sh.ds[M] * T)
	outc[:, M] .= tmp[2] .+ sh.flow.F(tmp[1]) .* sh.ds[M] * dT .- dxc[:, 1]

	# add constraint
	out[end] = (sh.section(x .+ δ .* dx) - sh.section(x)) / δ

	return out
end

# shooting functional, this allows for Array state space
function (sh::ShootingProblem)(x::BorderedArray)
	# period of the cycle
	T = extractPeriodShooting(x)
	M = length(x.u)

	# variable to hold the computed result
	out = BorderedArray{typeof(x.u), typeof(x.p)}(similar(x.u), typeof(x.p)(0))

	for ii in 1:M-1
		# yi = Flow(xi, dsi⋅T) - x[i+1]
		out.u[ii] .= sh.flow(x.u[ii], sh.ds[ii] * T) .- x.u[ii+1]
	end
	# ym = Flow(xm, sm T) - x1
	out.u[M] .= sh.flow(x.u[M], sh.ds[M] * T) .- x.u[1]

	# add constraint
	out.p = sh.section(x)

	return out
end

# jacobian of the shooting functional, this allows for Array state space
function (sh::ShootingProblem)(x::BorderedArray, dx::BorderedArray; δ = 1e-8)
	# period of the cycle
	dT = extractPeriodShooting(dx)
	T  = extractPeriodShooting(x)
	M = length(x.u)

	# variable to hold the computed result
	out = BorderedArray{typeof(x.u), typeof(x.p)}(similar(x.u), typeof(x.p)(0))

	for ii = 1:M-1
		# call jacobian of the flow
		tmp = sh.flow(x.u[ii], dx.u[ii], sh.ds[ii] * T)
		out.u[ii] .= tmp[2] .+ sh.flow.F(tmp[1]) .* sh.ds[ii] * dT .- dx.u[ii+1]
	end

	# ym = Flow(xm, sm T) - x1
	tmp = sh.flow(x.u[M], dx.u[M], sh.ds[M] * T)
	out.u[M] .= tmp[2] .+ sh.flow.F(tmp[1]) .* sh.ds[M] * dT .- dx.u[1]

	# add constraint
	x_tmp = similar(x.u); copyto!(x_tmp, x.u)
	axpy!(δ , dx.u, x_tmp)
	out.p = (sh.section(BorderedArray(x_tmp, T + δ * dT)) - sh.section(x)) / δ

	return out
end

function _getMax(prob::ShootingProblem, x::AbstractVector; ratio = 1)
	# this function extracts the amplitude of the cycle
	T = extractPeriodShooting(x)
	M = length(prob.ds)
	N = div(length(x) - 1, M)
	xv = @view x[1:end-1]
	xc = reshape(xv, N, M)

	# we could use @views but then Sundials might complain
	sol = prob.flow.flowFull((xc[:, 1]), T)
	mx = @views maximum(sol[1:div(N, ratio), :], dims = 1)
end

function getAmplitude(prob::ShootingProblem, x::AbstractVector; ratio = 1)
	mx = _getMax(prob, x; ratio = ratio)
	return maximum(mx) - minimum(mx)
end

function getMaximum(prob::ShootingProblem, x::AbstractVector; ratio = 1)
	mx = _getMax(prob, x; ratio = ratio)
	return maximum(mx)
end

function sectionShooting(x::AbstractVector, po::AbstractMatrix, p, F)
	res = 1.0
	N, M = size(po)
	xv = @view x[1:end-1]
	xc = reshape(xv, N, M)

	for ii = 1:size(x, 2)
		res *= @views dot(xc[:, ii], F(po[:, ii], p)) .- dot(po[:, ii], F(po[:, ii], p))
	end
	res
end

####################################################################################################
"""
pb = PoincareShootingProblem(flow::Flow, M, section)

This composite type implements the Poincaré Shooting method to locate periodic orbits, basically using Poincaré return maps. The arguments are as follows
- `flow::Flow`: implements the flow of the Cauchy problem though the structure `Flow`.
- `M`: the number of return maps. If `M==1`, then the simple shooting is implemented and the multiple one otherwise.
- `section`: implements a Poincaré section condition. The evaluation `section(x)` must return a scalar number where `x` is a guess for the periodic orbit when `M=1`. Otherwise, one must implement a function `section(out, x)` which populates `out` with the `M` hyperplanes.

## Simplified constructors
A simpler way is to create a functional is `pb = PoincareShootingProblem(F, p, prob::ODEProblem, alg, section; kwargs...)` for simple shooting or `pb = PoincareShootingProblem(F, p, prob::ODEProblem, alg, M::Int, section; kwargs...)` for multiple shooting . Here `F` is the vector field, `p` is a parameter (to be passed to the vector and the flow), `prob` is an `ODEProblem` which is used to create a flow using the ODE solver `alg` (for example `Tsit5()`). Finally, the arguments `kwargs` are passed to the ODE solver defining the flow. Look at `DifferentialEquations.jl` for more information.

Another convenient call is `pb = PoincareShootingProblem(F, p, prob::ODEProblem, alg, normals::AbstractVector, centers::AbstractVector; kwargs...)` where `normals` (resp. `centers`) is a list of normals (resp. centers) which define a list of hyperplanes `\\Sigma_i`. These hyperplanes are used to define partial Poincaré return maps. See docs for more information.

## Computing the functionals
You can then call `pb(orbitguess)` to apply the functional to a guess. Note that `orbitguess` must be of size M * N where N is the number of unknowns in the state space and `M` is the number of Poincaré maps. Another accepted `guess` is such that `guess[i]` is the state of the orbit on the `i`th section. This last form allows for non-vector state space which can be convenient for 2d problems for example.

A functional, hereby called `G` encodes this shooting problem. For example, the following methods are available:
- `pb(orbitguess)` evaluates the functional G on `orbitguess`
- `pb(orbitguess, du)` evaluates the jacobian `dG(orbitguess).du` functional at `orbitguess` on `du`

!!! tip "Tip"
    You can use the function `getPeriod(pb::PoincareShootingProblem, sol)` to get the period of the solution `sol`
"""
@with_kw struct PoincareShootingProblem{Tf <: Flow, Tsection} <: AbstractShootingProblem
	flow::Tf					# should be a Flow{TF, Tf, Td}
	M::Int64 = 1				# number of Poincaré sections
	section::Tsection			# Poincaré sections
end

# simple Poincaré shooting
function PoincareShootingProblem(F, prob::ODEProblem, alg, section; interp_points = 50, kwargs...)
	pSection(u, t, integrator) = section(u) .* (integrator.iter > 1)
	# we put the nothing option to have an upcrossing
	cb = ContinuousCallback(pSection, terminate!, nothing, interp_points = interp_points)
	return PoincareShootingProblem(Flow(F, prob, alg, callback=cb; kwargs...), 1, section)
end

# this function takes into account a parameter passed to the vector field
# simple shooting
function PoincareShootingProblem(F, p, prob::ODEProblem, alg, section; interp_points = 50, kwargs...)
	pSection(u, t, integrator) = section(u) * (integrator.iter > 1)
	# we put nothing option to have an upcrossing
	cb = ContinuousCallback(pSection, terminate!, nothing, interp_points = interp_points)
	return PoincareShootingProblem(Flow(F, p, prob, alg; callback = cb, kwargs...), 1, section)
end

# this function takes into account a parameter passed to the vector field
# multiple shooting
function PoincareShootingProblem(F, p, prob::ODEProblem, alg, M, section; interp_points = 50, kwargs...)
	if M==1
		return PoincareShootingProblem(F, p, prob::ODEProblem, alg, section; kwargs...)
	end

	pSection(out, u, t, integrator) = section(out, u) .* (integrator.iter > 1)
	affect!(integrator, idx) = terminate!(integrator)
	# we put the nothing option to have an upcrossing
	cb = VectorContinuousCallback(pSection, affect!, M;interp_points = interp_points, affect_neg! = nothing)
	return PoincareShootingProblem(Flow(F, p, prob, alg; callback = cb, kwargs...), M, section)
end

function getPeriod(prob::PoincareShootingProblem, x)
	# this function extracts the period of the cycle
	@warn "Only done for M=1"
	prob.flow.flowFull(x, Inf64).t[end]
end

# Poincaré shooting functional
function (sh::PoincareShootingProblem)(x::AbstractVector)
	# period of the cycle
	M = sh.M
	N = div(length(x), M)

	# reshape the period orbit guess
	xc = reshape(x, N, M)

	# variable to hold the computed result
	out = similar(x)
	outc = reshape(out, N, M)

	for ii in 1:M-1
		# yi = Flow(xi) - x[i+1]
		# here we compute the flow up to an infinite time hoping we intersect a Poincare section before :D
		@views outc[:, ii] .= sh.flow(xc[:, ii], Inf64) .- xc[:, ii+1]
	end
	# ym = Flow(xm) - x1
	@views outc[:, M] .= sh.flow(xc[:, M], Inf64) .- xc[:, 1]

	return out
end

# jacobian of the shooting functional
function (sh::PoincareShootingProblem)(x::AbstractVector, dx::AbstractVector; δ = 1e-8)
	return (sh(x .+ δ .* dx) .- sh(x)) ./ δ
end
####################################################################################################
# Poincare shooting based on Sánchez, J., M. Net, B. Garcı́a-Archilla, and C. Simó. “Newton–Krylov Continuation of Periodic Orbits for Navier–Stokes Flows.” Journal of Computational Physics 201, no. 1 (November 20, 2004): 13–33. https://doi.org/10.1016/j.jcp.2004.04.018.

function sectionHyp(out, x, normals, centers)
	for ii = 1:length(normals)
		out[ii] = dot(normals[ii], x .- centers[ii])
	end
	out
end

# this composite type encodes a set of hyperplanes which are used as Poincaré sections
struct HyperplaneSections
	M				# number of hyperplanes
	normals 		# normals to define hyperplanes
	centers 		# representative point on each hyperplane
	indices 		# indices to be removed in the operator Ek

	function HyperplaneSections(normals, centers)
		M = length(normals)
		indices = zeros(Int64, M)
		for ii=1:M
			indices[ii] = argmax(abs.(normals[ii]))
		end
		return new(M, normals, centers, indices)
	end
end

function (hyp::HyperplaneSections)(out, u)
	sectionHyp(out, u, hyp.normals, hyp.centers)
end

# Operateur Rk from the paper above
function R!(hyp::HyperplaneSections, out, x::AbstractVector, ii::Int)
	k = hyp.indices[ii]
	@views out[1:k-1] .= x[1:k-1]
	@views out[k:end] .= x[k+1:end]
	return out
end

function R(hyp::HyperplaneSections, x::AbstractVector, ii::Int)
	out = similar(x, length(x) - 1)
	R!(hyp, out, x, ii)
end

# Operateur Ek from the paper above
function E!(hyp::HyperplaneSections, out, xbar::AbstractVector, ii::Int)
	@assert length(xbar) == length(hyp.normals[1]) - 1 "Wrong size for the projector / expansion operators, length(xbar) = $(length(xbar)) and length(normal) = $(length(hyp.normals[1]))"
	k = hyp.indices[ii]
	nbar = R(hyp, hyp.normals[ii], ii)
	xcbar = R(hyp, hyp.centers[ii], ii)
	coord_k = hyp.centers[ii][k] - dot(nbar, xbar .- xcbar) / hyp.normals[ii][k]

	@views out[1:k-1] .= xbar[1:k-1]
	@views out[k+1:end] .= xbar[k:end]
	out[k] = coord_k

	return out
end

function E(hyp::HyperplaneSections, xbar::AbstractVector, ii::Int)
	out = similar(xbar, length(xbar) + 1)
	E!(hyp, out, xbar, ii)
end

struct HyperplanePoincareShootingProblem <: AbstractShootingProblem
	psh::PoincareShootingProblem
end

getPeriod(hpsh::HyperplanePoincareShootingProblem, x) = getPeriod(hpsh.psh)

function PoincareShootingProblem(F, p, prob::ODEProblem, alg, normals::AbstractVector, centers::AbstractVector; interp_points = 50, kwargs...)
	hyp = HyperplaneSections(normals, centers)
	pSection(out, u, t, integrator) = hyp(out, u) * (integrator.iter > 1)
	affect!(integrator, idx) = terminate!(integrator)
	# we put nothing option to have an upcrossing
	cb = VectorContinuousCallback(pSection, affect!, hyp.M; interp_points = interp_points, affect_neg! = nothing)
	return HyperplanePoincareShootingProblem(PoincareShootingProblem(Flow(F, p, prob, alg; callback = cb, kwargs...), hyp.M, hyp))
end

# Poincaré (multiple) shooting with hyperplanes parametrization
function (hpsh::HyperplanePoincareShootingProblem)(x_bar::AbstractVector; verbose = false)
	sh = hpsh.psh
	M = sh.M
	Nm1 = div(length(x_bar), M)

	# reshape the period orbit guess
	x_barc = reshape(x_bar, Nm1, M)

	# variable to hold the computed result
	xc = similar(x_bar, Nm1 + 1, M)
	outc = similar(xc)

	# we extend the state space to be able to call the flow, so we fill xc
	#TODO create the projections on the fly
	for ii=1:M
		E!(sh.section, view(xc, :, ii), view(x_barc, :, ii), ii)
	end

	for ii in 1:M-1
		# yi = Flow(xi) - x[i+1]
		@views outc[:, ii] .= xc[:, ii+1].- sh.flow(xc[:, ii], Inf64)
	end
	# ym = Flow(xm) - x1
	@views outc[:, M] .= xc[:, 1] .- sh.flow(xc[:, M], Inf64)

	# build the array to be returned
	out_bar = similar(x_bar)
	out_barc = reshape(out_bar, Nm1, M)
	for ii=1:M-1
		R!(hpsh.psh.section, view(out_barc, :, ii), view(outc, :, ii), ii+1)
	end
	R!(hpsh.psh.section, view(out_barc, :, M), view(outc, :, M), 1)

	return out_bar
end

# jacobian of the shooting functional
function (sh::HyperplanePoincareShootingProblem)(x::AbstractVector, dx::AbstractVector; δ = 1e-8)
	return (sh(x .+ δ .* dx) .- sh(x)) ./ δ
end

####################################################################################################
# if we use the same code as for newton (see below) in continuation, it is difficult to tell the eigensolver not to use the jacobian but instead the monodromy matrix. So we have to use a dedicated composite type for the jacobian to handle this case.

struct ShootingJacobian{Tpb <: AbstractShootingProblem, Torbitguess}
	pb::Tpb
	x::Torbitguess
end

(shjac::ShootingJacobian)(dx) = shjac.pb(shjac.x, dx)

####################################################################################################
# newton wrapper
"""
	newton(prob::T, orbitguess, options::NewtonPar; kwargs...) where {T <: AbstractShootingProblem}

This is the Newton Solver for computing a periodic orbit using Shooting method.

# Output:
- solution
- history of residuals
- flag of convergence
- number of iterations
"""
function newton(prob::T, orbitguess, options::NewtonPar; kwargs...) where {T <: AbstractShootingProblem}
	return newton(x -> prob(x),
			x -> (dx -> prob(x, dx)),
			orbitguess,
			options; kwargs...)
end

"""
	newton(prob::T, orbitguess, options::NewtonPar, defOp::DeflationOperator; kwargs...) where {T <: AbstractShootingProblem}

This is the deflated Newton Solver for computing a periodic orbit using Shooting method.

# Output:
- solution
- history of residuals
- flag of convergence
- number of iterations
"""
function newton(prob::Tpb, orbitguess, options::NewtonPar, defOp::DeflationOperator{T, Tf, vectype}; kwargs...) where {Tpb <: AbstractShootingProblem, T, Tf, vectype}
	return newton(x -> prob(x),
			x -> (dx -> prob(x, dx)),
			orbitguess,
			options, defOp; kwargs...)
end

####################################################################################################
# Continuation

function continuationPOShooting(prob, orbitguess, p0::Real, _contParams::ContinuationPar, linearAlgo::AbstractBorderedLinearSolver; printPeriod = true, kwargs...)

	contParams = check(_contParams)

	options = contParams.newtonOptions

	if contParams.computeEigenValues
		contParams = @set contParams.newtonOptions.eigsolver = FloquetQaDShooting(contParams.newtonOptions.eigsolver)
	end

	if (prob(p0) isa PoincareShootingProblem) || (prob(p0) isa HyperplanePoincareShootingProblem)
		if printPeriod
			printSolutionPS = (x, p) -> getPeriod(prob(p), x)
			return continuation(
				(x, p) -> prob(p)(x),
				# (x, p) -> (dx -> prob(p)(x, dx)),
				(x, p) -> ShootingJacobian(prob(p), x),
				orbitguess, p0,
				contParams, linearAlgo;
				printSolution = printSolutionPS,
				kwargs...)
		end
	end

	return continuation(
		(x, p) -> prob(p)(x),
		(x, p) -> ShootingJacobian(prob(p), x),
		orbitguess, p0,
		contParams, linearAlgo;
		kwargs...)
end

"""
	continuationPOShooting(prob, orbitguess, p0::Real, contParams::ContinuationPar; printPeriod = true, kwargs...)

This is the continuation routine for computing a periodic orbit using a functional G based on a Shooting method.

# Arguments
- `p -> prob(p)` is a function or family such that `prob(p)::AbstractShootingProblem` encodes the functional G
- `orbitguess` a guess for the periodic orbit. For the type of `orbitguess`, please see the information concerning `ShootingProblem` and `PoincareShootingProblem`.
- `p0` initial parameter, must be a real number
- `contParams` same as for the regular `continuation` method
- `printPeriod` in the case of Poincaré Shooting, plot the period of the cycle.
"""
function continuationPOShooting(prob, orbitguess, p0::Real, contParams::ContinuationPar; linearAlgo = BorderingBLS(), printPeriod = true, kwargs...)
	_linearAlgo = @set linearAlgo.solver = contParams.newtonOptions.linsolver
	return continuationPOShooting(prob, orbitguess, p0, contParams, _linearAlgo; printPeriod = printPeriod, kwargs...)
end