text:
"""
Driver to run TRON on GPU. This should be called from a kernel.
"""
@inline function tron_kernel(n::Int, shift::Int, max_feval::Int, max_minor::Int, gtol::Float64, scale::Float64, use_polar::Bool,
                     x, xl,
                     xu,
                     param,
                     YffR::Float64, YffI::Float64,
                     YftR::Float64, YftI::Float64,
                     YttR::Float64, YttI::Float64,
                     YtfR::Float64, YtfI::Float64,
                     I, J
                     )

    tx = J
    g = @localmem Float64 (4,)
    xc = @localmem Float64 (4,)
    s = @localmem Float64 (4,)
    wa = @localmem Float64 (4,)
    wa1 = @localmem Float64 (4,)
    wa2 = @localmem Float64 (4,)
    wa3 = @localmem Float64 (4,)
    wa4 = @localmem Float64 (4,)
    wa5 = @localmem Float64 (4,)
    wa = @localmem Float64 (4,)
    gfree = @localmem Float64 (4,)
    dsave = @localmem Float64 (4,)
    indfree = @localmem Int (4,)
    iwa = @localmem Int (4,)
    isave = @localmem Int (4,)

    A = @localmem Float64 (4,4)
    B = @localmem Float64 (4,4)
    L = @localmem Float64 (4,4)

    if tx <= n
        @inbounds for j=1:n
            A[tx,j] = 0.0
            B[tx,j] = 0.0
            L[tx,j] = 0.0
        end
    end
    @synchronize

    task = 0
    status = 0

    delta = 0.0
    fatol = 0.0
    frtol = 1e-12
    fmin = -1e32
    cgtol = 0.1
    cg_itermax = n

    f = 0.0
    nfev = 0
    ngev = 0
    nhev = 0
    minor_iter = 0
    search = true

    while search

        # [0|1]: Evaluate function.

        if task == 0 || task == 1
            if use_polar
                f = eval_f_polar_kernel(n, shift, scale, x, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI, I, J)
            else
                f = eval_f_kernel(n, scale, x, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI, I, J)
            end
            nfev += 1
            if nfev >= max_feval
                search = false
            end

        end

        # [2] G or H: Evaluate gradient and Hessian.

        if task == 0 || task == 2
            if use_polar
                eval_grad_f_polar_kernel(n, shift, scale, x, g, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI, I, J)
                eval_h_polar_kernel(n, shift, scale, x, A, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI, I, J)
            else
                eval_grad_f_kernel(n, scale, x, g, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI, I, J)
                eval_h_kernel(n, scale, x, A, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI, I, J)
            end
            ngev += 1
            nhev += 1
            minor_iter += 1
        end

        # Initialize the trust region bound.

        if task == 0
            gnorm0 = dnrm2(n, g, 1, J)
            delta = gnorm0
        end

        # Call Tron.

        if search
            delta, task = ExaTronKernels.dtron(n, x, xl, xu, f, g, A, frtol, fatol, fmin, cgtol,
                                cg_itermax, delta, task, B, L, xc, s, indfree, gfree,
                                isave, dsave, wa, iwa, wa1, wa2, wa3, wa4, wa5, J)
        end

        # [3] NEWX: a new point was computed.

        if task == 3
            gnorm_inf = ExaTronKernels.dgpnorm(n, x, xl, xu, g, J)

            if gnorm_inf <= gtol
                task = 4
            end

            if minor_iter >= max_minor
                status = 1
                search = false
            end
        end

        # [4] CONV: convergence was achieved.
        # [10] : warning fval is less than fmin

        if task == 4 || task == 10
            search = false
        end
    end

    @synchronize

    return status, minor_iter
end


text:
using LinearAlgebra
using SparseArrays
using Logging
BLAS.set_num_threads(1)  # 1 thread only

using QPSReader

import Tulip
using UnitBlockAngular

function run_rmp_tulip_uba(frmp)
    finst = split(frmp, '/')[end]

    # Check if DER or TSSP instance
    if finst[1:3] == "DER"
        R = parse(Int, split(finst, '_')[3])
    else
        # TSSP
        R = parse(Int, split(finst, r"_|\.")[2])
    end

    # Read RMP
    qps = with_logger(Logging.NullLogger()) do
        readqps(frmp)
    end

    # Problem stats:
    #   * linking constraints
    #   * linking variables
    #   * Number of sub-problems
    #   * Density of lower blocks
    M = qps.ncon
    N = qps.nvar
    m0 = M - R

    nz = 0
    blockidx = zeros(Int, N)
    for (i, j, v) in zip(qps.arows, qps.acols, qps.avals)
        if 1 <= i <= m0
            nz += 1
        else
            @assert blockidx[j] == 0
            blockidx[j] = i - m0
        end
    end
    n0 = sum(blockidx .== 0)

    # Sort columns
    perm = sortperm(blockidx)
    iperm = invperm(perm)

    # Re-order columns and load model
    acols_ = copy(qps.acols)
    for (k, j) in enumerate(qps.acols)
        acols_[k] = iperm[j]
    end

    tlp = Tulip.Model{Float64}()
    Tulip.load_problem!(tlp.pbdata,
        qps.name,
        true, qps.c[perm], qps.c0,
        sparse(qps.arows, acols_, qps.avals, qps.ncon, qps.nvar),
        qps.lcon, qps.ucon,
        qps.lvar[perm], qps.uvar[perm],
        qps.connames, qps.varnames[perm]
    )

    # Set parameters
    tlp.params.OutputLevel = 1
    tlp.params.Presolve = 0

    # Set custom KKT options
    tlp.params.MatrixOptions = Tulip.TLA.MatrixOptions(UnitBlockAngularMatrix, m0=m0, n0=n0, n=N-n0, R=R)
    tlp.params.KKTOptions = Tulip.KKT.SolverOptions(UnitBlockAngularFactor)

    # Solve model
    t = @elapsed Tulip.optimize!(tlp)

    println("Solve time: $t")
    println("Barrier iterations: ", Tulip.get_attribute(tlp, Tulip.BarrierIterations()))

    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    run_rmp_tulip_uba(ARGS[1])
end

ext: jl

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
# Unit to calculate solar geometry including solar elevation, declination,
#  azimuth etc using TSolar class. Data are hidden. 03/15/00 SK
# - 1st Revision 10/10/00: Changed to dealing only upto the top of the canopy. Radiation transter within the canopy is now a separate module.
# - added functins to calculate global radiation, atmospheric emissivity, etc as in Spitters et al. (1986), 3/18/01, SK
# 24Dec03, SK
# - added a function to calculate day length based on Campbell and Norman (1998) p 170,
# - added overloads to SetVal
# 2Aug04, SK
# - Translated to C++ from Delphi
# - revised some functions according to "An introduction to solar radiaiton" by Iqbal (1983)
# - (To Do) Add algorithms for instantaneous diffuse and direct radiation predictions from daily global solar radiation for given time
# - (To Do) This can be done by first applying sinusoidal model to the daily data to simulate hourly global solar radiation
# - (To Do) Then the division model of diffuse and direct radiations was applied
# - added direct and diffuse separation functions according to Weiss and Norman (1985), 3/16/05

@system Location begin
    lat: latitude => 36u"°" ~ preserve(u"°", parameter)
    long: longitude => 128u"°" ~ preserve(u"°", parameter)
    alt: altitude => 20u"m" ~ preserve(u"m", parameter)
end

import Dates

@system Sun begin
    #TODO override Location
    loc(context): location ~ ::Location #(override)
    calendar ~ ::Calendar(override)
    weather ~ ::Weather(override)

    # @derive time? -- takes account different Julian day conventions (03-01 vs. 01-01)
    t(calendar.time): datetime ~ track::datetime
    d(t): day => Dates.dayofyear(t) ~ track::int(u"d")
    h(t): hour => Dates.hour(t) ~ track::int(u"hr")

    ϕ(loc.lat): latitude ~ track(u"°") # DO NOT convert to radians for consistency
    λ(loc.long): longitude ~ track(u"°") # leave it as in degrees, used only once for solar noon calculation
    alt(loc.alt): altitude ~ track(u"m")

    solrad(weather.solrad): solar_radiation ~ track(u"W/m^2")
    τ: transmissivity => 0.5 ~ preserve(parameter) # atmospheric transmissivity, Goudriaan and van Laar (1994) p 30

    #####################
    # Solar Coordinates #
    #####################

    #HACK always use degrees for consistency and easy tracing
    #FIXME pascal version of LightEnv uses iqbal()
    δ(declination_angle_spencer): declination_angle ~ track(u"°")

    # Goudriaan 1977
    declination_angle_goudriaan(d) => begin
        g = 2pi * (d + 10u"d") / 365u"d"
        -23.45u"°" * cos(g)
    end ~ track(u"°")

    # Resenberg, blad, verma 1982
    declination_angle_resenberg(d) => begin
        g = 2pi * (d - 172u"d") / 365u"d"
        23.5u"°" * cos(g)
    end ~ track(u"°")

    # Iqbal (1983) Pg 10 Eqn 1.3.3, and sundesign.com
    declination_angle_iqbal(d) => begin
        g = 2pi * (d + 284u"d") / 365u"d"
        23.45u"°" * sin(g)
    end ~ track(u"°")

    # Campbell and Norman, p168
    declination_angle_campbell(d) => begin
        a = deg2rad(356.6 + 0.9856u"d^-1" * d)
        b = deg2rad(278.97 + 0.9856u"d^-1" * d + 1.9165sin(a))
        asind(0.39785sin(b))
    end ~ track(u"°")

    # Spencer equation, Iqbal (1983) Pg 7 Eqn 1.3.1. Most accurate among all
    declination_angle_spencer(d) => begin
        # gamma: day angle
        g = 2pi * (d - 1u"d") / 365u"d"
        0.006918 - 0.399912cos(g) + 0.070257sin(g) - 0.006758cos(2g) + 0.000907sin(2g) -0.002697cos(3g) + 0.00148sin(3g)
    end ~ track(u"rad")

    dph: degree_per_hour => 360u"°" / 24u"hr" ~ preserve(u"°/hr")

    # LC is longitude correction for Light noon, Wohlfart et al, 2000; Campbell & Norman 1998
    LC(λ, dph): longitude_correction => begin
        # standard meridian for pacific time zone is 120 W, Eastern Time zone : 75W
        # LC is positive if local meridian is east of standard meridian, i.e., 76E is east of 75E
        #standard_meridian = -120
        meridian = round(u"hr", λ / dph) * dph
        #FIXME use standard longitude sign convention
        #(long - meridian) / dph
        #HACK this assumes inverted longitude sign that MAIZSIM uses
        (meridian - λ) / dph
    end ~ track(u"hr")

    ET(d): equation_of_time => begin
        f = (279.575 + 0.9856u"d^-1" * d)*u"°"
        (-104.7sin(f) + 596.2sin(2f) + 4.3sin(3f) - 12.7sin(4f) -429.3cos(f) - 2.0cos(2f) + 19.3cos(3f)) / (60 * 60)
    end ~ track(u"hr")

    solar_noon(LC, ET) => 12u"hr" - LC - ET ~ track(u"hr")

    # θs: zenith angle
    hour_angle_at(ϕ, δ; θs(u"°")) => begin
        # this value should never become negative because -90 <= latitude <= 90 and -23.45 < decl < 23.45
        #HACK is this really needed for crop models?
        # preventing division by zero for N and S poles
        #denom = fmax(denom, 0.0001)
        # sunrise/sunset hour angle
        #TODO need to deal with lat_bound to prevent tan(90)?
        #lat_bound = radians(68)? radians(85)?
        # cos(h0) at cos(theta_s) = 0 (solar zenith angle = 90 deg == elevation angle = 0 deg)
        #-tan(ϕ) * tan(δ)
        c = (cos(θs) - sin(ϕ) * sin(δ)) / (cos(ϕ) * cos(δ))
        # c > 1: in the polar region during the winter, sun does not rise
        # c < -1: white nights during the summer in the polar region
        c = clamp(c, -1, 1)
        acosd(c)
    end ~ call(u"°")

    hour_angle_at_horizon(hour_angle_at) => hour_angle_at(90u"°") ~ track(u"°")

    # from Iqbal (1983) p 16
    half_day_length(hour_angle_at_horizon, dph) => (hour_angle_at_horizon / dph) ~ track(u"hr")
    day_length(half_day_length) => 2half_day_length ~ track(u"hr")

    sunrise(solar_noon, half_day_length) => (solar_noon - half_day_length) ~ track(u"hr")
    sunset(solar_noon, half_day_length) => (solar_noon + half_day_length) ~ track(u"hr")

    hour_angle(h, solar_noon, dph) => ((h - solar_noon) * dph) ~ track(u"°")

    αs(h=hour_angle, δ, ϕ): elevation_angle => begin
        #FIXME When time gets the same as solarnoon, this function fails. 3/11/01 ??
        asind(cos(h) * cos(δ) * cos(ϕ) + sin(δ) * sin(ϕ))
    end ~ track(u"°")

    ts(αs): positive_elevation_angle ~ track(u"°", min=0)

    θs(αs): zenith_angle => (90u"°" - αs) ~ track(u"°")

    # The solar azimuth angle is the angular distance between due South and the
    # projection of the line of sight to the sun on the ground.
    # View point from south, morning: +, afternoon: -
    # See An introduction to solar radiation by Iqbal (1983) p 15-16
    # Also see https://www.susdesign.com/sunangle/
    ϕs(αs, δ, ϕ): azimuth_angle => begin
        acosd((sin(δ) - sin(αs) * sin(ϕ)) / (cos(αs) * cos(ϕ)))
    end ~ track(u"°")

    ###################
    # Solar Radiation #
    ###################

    # atmospheric pressure in kPa (default = 100 kPa?)
    p(altitude): atmospheric_pressure => begin
        # campbell and Norman (1998), p 41
        101.3exp(-altitude / 8200u"m")
    end ~ track(u"kPa")

    m(p, ts): optical_air_mass_number => begin
        #FIXME check 101.3 is indeed in kPa
        #iszero(t_s) ? 0. : p / (101.3u"kPa" * sin(t_s))
        p / (101.3u"kPa" * sin(ts))
    end ~ track

    SC: solar_constant => 1370 ~ preserve(u"W/m^2", parameter)

    # Campbell and Norman's global solar radiation, this approach is used here
    insolation(ts, d, SC) => begin
        # solar constant, Iqbal (1983)
        #FIXME better to be 1361 or 1362 W/m-2?
        g = 2pi * (d - 10u"d") / 365u"d"
        SC * sin(ts) * (1 + 0.033cos(g))
    end ~ track(u"W/m^2")

    directional_solar_radiation(Fdir, solar_radiation) => begin
        Fdir * solar_radiation
    end ~ track(u"W/m^2")

    diffusive_solar_radiation(Fdif, solar_radiation) => begin
        Fdif * solar_radiation
    end ~ track(u"W/m^2")

    Fdir(τ, m): directional_coeff => begin
        # Goudriaan and van Laar's global solar radiation
        #FIXME should be goudriaan() version
        goudriaan(τ) = τ * (1 - diffusive_coeff)
        #FIXME: check if equation is same as campbell()
        # Takakura (1993), p 5.11
        takakura(τ, m) = τ^m
        # Campbell and Norman (1998), p 173
        campbell(τ, m) = τ^m
        campbell(τ, m)
    end ~ track

    # Fdif: Fraction of diffused light
    Fdif(τ, m): diffusive_coeff => begin
        # Goudriaan and van Laar's global solar radiation
        goudriaan(τ) = begin
            # clear sky : 20% diffuse
            if τ >= 0.7
                0.2
            # cloudy sky: 100% diffuse
            elseif τ <= 0.3
                1
            # inbetween
            else
                1.6 - 2τ
            end
        end
        # Takakura (1993), p 5.11
        takakura(τ, m) = (1 - τ^m) / (1 - 1.4log(τ)) / 2
        # Campbell and Norman (1998), p 173
        campbell(τ, m) = 0.3(1 - τ^m)
        campbell(τ, m)
    end ~ track

    directional_fraction(Fdif, Fdir) => (1 / (1 + Fdif/Fdir)) ~ track
    diffusive_fraction(Fdir, Fdif) => (1 / (1 + Fdir/Fdif)) ~ track

    # PARfr
    #TODO better naming: extinction? transmitted_fraction?
    PARfr(τ): photosynthetic_coeff => begin
        #if self.elevation_angle <= 0:
        #    0
        #TODO: implement Weiss and Norman (1985), 3/16/05
        weiss() = nothing
        # Goudriaan and van Laar (1994)
        goudriaan(τ) = begin
            # clear sky (τ >= 0.7): 45% is PAR
            if τ >= 0.7
                0.45
            # cloudy sky (τ <= 0.3): 55% is PAR
            elseif τ <= 0.3
                0.55
            else
                0.625 - 0.25τ
            end
        end
        goudriaan(τ)
    end ~ track

    Q: photosynthetic_active_radiation_conversion_factor => begin
        # 4.55 is a conversion factor from W to photons for solar radiation, Goudriaan and van Laar (1994)
        # some use 4.6 i.e., Amthor 1994, McCree 1981, Challa 1995.
        4.6
    end ~ preserve(u"μmol/J", parameter)

    # PARtot: total PAR (umol m-2 s-1) on horizontal surface (PFD)
    PARtot(solar_radiation, PARfr, Q): photosynthetic_active_radiation_total => begin
        # conversion factor from W/m2 to PFD (umol m-2 s-1) for PAR waveband (median 550 nm of 400-700 nm) of solar radiation,
        # see Campbell and Norman (1994) p 149
        solar_radiation * PARfr * Q
    end ~ track(u"μmol/m^2/s") # Quanta

    # PARdir
    PARdir(PARtot, directional_fraction): directional_photosynthetic_radiation => (directional_fraction * PARtot) ~ track(u"μmol/m^2/s") # Quanta

    # PARdif
    PARdif(PARtot, diffusive_fraction): diffusive_photosynthetic_radiation => (diffusive_fraction * PARtot) ~ track(u"μmol/m^2/s") # Quanta
end

@system SunController(Controller) begin
    calendar(context) ~ ::Calendar
    weather(context, calendar) ~ ::Weather
    s(context, calendar, weather): sun ~ ::Sun
end

plot_sun(v, d=3u"d"; kw...) = begin
    o = (
        :Calendar => (:init => ZonedDateTime(2007, 9, 1, tz"Asia/Seoul")),
        :Weather => (:store => loadwea(datapath("2007.wea"), tz"Asia/Seoul")),
    )
    r = simulate(SunController, stop=d, config=o, base=:sun)
    visualize(r, :time, v; kw...)
end

test_sun(d=3u"d") = foreach(v -> display(plot_sun(v, d)), [
    :δ, # declination_angle
    :αs, # elevation_angle
    :Fdir, # directional_coeff
    :Fdif, # diffusive_coeff
])


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id: condenser.mso 555 2008-07-18 19:01:13Z rafael $
#*--------------------------------------------------------------------
type condenser
	condenser()=begin
		PP=outers.PP
		NComp=outers.NComp
		[
			:(LI = Levelpercent_Initial),
			:(OutletLiquid.T = Initial_Temperature),
			:(OutletLiquid.z(1:NComp-1) = Initial_Composition(1:NComp-1)/sum(Initial_Composition)),
		],
		[
			"Initial level Percent","Initial Temperature","Initial Composition",
		],
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Components"
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight",
				:Hidden=>true
			)),(NComp)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"Low Flow",
				:Default=>1E-6,
				:Hidden=>true
			)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"No Flow",
				:Default=>0,
				:Hidden=>true
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Vapour Flow",
				:Valid=>["on", "off"],
				:Default=>"on",
				:Hidden=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"K factor for pressure drop",
				:Lower=>1E-8,
				:Default=>1E-3
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Initial liquid height in Percent",
				:Default=>0.70
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Initial Temperature of Condenser"
			)),
			fill(positive(Dict{Symbol,Any}(
				:Brief=>"Initial Liquid Composition",
				:Lower=>1E-6
			)),(NComp)),
			VesselVolume(Dict{Symbol,Any}(
				:Brief=>"Vessel Geometry",
				:Symbol=>" "
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Vapour inlet stream",
				:PosX=>0.13,
				:PosY=>0,
				:Symbol=>"_{in}^{Vapour}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream",
				:PosX=>0.35,
				:PosY=>1,
				:Symbol=>"_{out}^{Liquid}"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour outlet stream",
				:PosX=>0.54,
				:PosY=>0,
				:Symbol=>"_{out}^{Vapour}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Heat supplied",
				:Protected=>true,
				:PosX=>1,
				:PosY=>0.08,
				:Symbol=>"Q_{in}"
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Temperature  Indicator of Condenser",
				:Protected=>true,
				:PosX=>0.33,
				:PosY=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Level  Indicator of Condenser",
				:Protected=>true,
				:PosX=>0.43,
				:PosY=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Pressure  Indicator of Condenser",
				:Protected=>true,
				:PosX=>0.25,
				:PosY=>0
			)),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Molar Holdup in the tray",
				:Protected=>true
			)),(NComp)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar liquid holdup",
				:Protected=>true
			)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar vapour holdup",
				:Protected=>true
			)),
			energy(Dict{Symbol,Any}(
				:Brief=>"Total Energy Holdup on tray",
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Volume",
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar volume",
				:Protected=>true
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Inlet Vapour Mass Density",
				:Hidden=>true,
				:Symbol=>"\\rho"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop",
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P",
				:Protected=>true
			)),
			[
				:(InletVapour.F*vV = Kfactor *sqrt(Pdrop/rho)*"m^2"),
				:(InletVapour.F = zero_flow),
				:(diff(M) = InletVapour.F*InletVapour.z - OutletLiquid.F*OutletLiquid.z- OutletVapour.F*OutletVapour.z),
				:(diff(E) = InletVapour.F*InletVapour.h - OutletLiquid.F*OutletLiquid.h- OutletVapour.F*OutletVapour.h + InletQ),
				:(M = ML*OutletLiquid.z + MV*OutletVapour.z),
				:(E = ML*OutletLiquid.h + MV*OutletVapour.h - OutletVapour.P*Geometry.Vtotal),
				:(sum(OutletLiquid.z)=1.0),
				:(sum(OutletLiquid.z)=sum(OutletVapour.z)),
				:(vL = PP.LiquidVolume(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)),
				:(vV = PP.VapourVolume(OutletVapour.T, OutletVapour.P, OutletVapour.z)),
				:(rho = PP.VapourDensity(InletVapour.T, InletVapour.P, InletVapour.z)),
				:(PP.LiquidFugacityCoefficient(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)*OutletLiquid.z = PP.VapourFugacityCoefficient(OutletVapour.T, OutletVapour.P, OutletVapour.z)*OutletVapour.z),
				:(OutletLiquid.T = OutletVapour.T),
				:(OutletVapour.P = OutletLiquid.P),
				:(OutletLiquid.P = InletVapour.P - Pdrop),
				:(Geometry.Vtotal = ML*vL + MV*vV),
				:(ML * vL = Geometry.Vfilled),
				:(TI * "K" = OutletLiquid.T),
				:(PI * "atm" = OutletLiquid.P),
				:(LI*Geometry.Vtotal= Geometry.Vfilled),
			],
			[
				"","","Component Molar Balance","Energy Balance","Molar Holdup","Energy Holdup","Mol fraction normalisation","Mol fraction Constraint","Liquid Volume","Vapour Volume","Inlet Vapour Density","Chemical Equilibrium","Thermal Equilibrium","Mechanical Equilibrium","Pressure Drop","Geometry Constraint","Liquid Level","Temperature indicator","Pressure indicator","Level indicator",
			],
			[:PP,:NComp,:Mw,:low_flow,:zero_flow,:VapourFlow,:Kfactor,:Levelpercent_Initial,:Initial_Temperature,:Initial_Composition,],
			[:Geometry,:InletVapour,:OutletLiquid,:OutletVapour,:InletQ,:TI,:LI,:PI,:M,:ML,:MV,:E,:vL,:vV,:rho,:Pdrop,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Mw::Array{molweight}
	low_flow::flow_mol
	zero_flow::flow_mol
	VapourFlow::DanaSwitcher
	Kfactor::positive
	Levelpercent_Initial::positive
	Initial_Temperature::temperature
	Initial_Composition::Array{positive}
	Geometry::VesselVolume
	InletVapour::stream
	OutletLiquid::liquid_stream
	OutletVapour::vapour_stream
	InletQ::power
	TI::control_signal
	LI::control_signal
	PI::control_signal
	M::Array{mol}
	ML::mol
	MV::mol
	E::energy
	vL::volume_mol
	vV::volume_mol
	rho::dens_mass
	Pdrop::press_delta
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export condenser
function set(in::condenser)
	Mw = PP.MolecularWeight()
	 low_flow = 1E-6 * "kmol/h"
	 zero_flow = 0 * "kmol/h"

end
function initial(in::condenser)
	addEquation(1)
	addEquation(2)
	addEquation(3)
end
function setEquationFlow(in::condenser)
	let switch=VapourFlow
		if VapourFlow==InletVapour.F < low_flow
			set(switch,"off")
		end
		if VapourFlow==InletVapour.P > OutletLiquid.P
			set(switch,"on")
		end
		if switch=="on"
			addEquation(1)
		elseif switch=="off"
			addEquation(2)
		end
	end
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
end
function atributes(in::condenser,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Condenser"
	fields[:Brief]="Model of a  dynamic condenser with control."
	fields[:Info]="== ASSUMPTIONS ==
* perfect mixing of both phases;
* thermodynamics equilibrium.

== SPECIFY ==
* the InletVapour stream;
* the outlet flows: OutletVapour.F and OutletLiquid.F;
* the InletQ (the model requires an energy stream, also you can use a controller for setting the heat duty using the heat_flow model).

== OPTIONAL ==
* the condenser model has three control ports
** TI OutletLiquid Temperature Indicator;
** PI OutletLiquid Pressure Indicator;
** LI Level Indicator of Condenser;

== INITIAL CONDITIONS ==
* Initial_Temperature :  the condenser temperature (OutletLiquid.T);
* Levelpercent_Initial : the condenser liquid level in percent (LI);
* Initial_Composition : (NoComps) OutletLiquid compositions.
"
	drive!(fields,_)
	return fields
end
condenser(_::Dict{Symbol,Any})=begin
	newModel=condenser()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
"""
### Kinetic matrix elements
"""
function kinetic( a::PGBF, b::PGBF )
    K = kinetic(a.expn, a.center, a.power, b.expn, b.center, b.power)
    return a.NORM*b.NORM*K
end

function kinetic( aexpn::Float64, acenter::Tuple3F64, apower::Tuple3I64,
                  bexpn::Float64, bcenter::Tuple3F64, bpower::Tuple3I64 )

    ax, ay, az = acenter[1],acenter[2],acenter[3]
    bx, by, bz = bcenter[1],bcenter[2],bcenter[3]

    aI,aJ,aK = apower[1],apower[2],apower[3]
    bI,bJ,bK = bpower[1],bpower[2],bpower[3]

    overlap0 = overlap( aexpn, acenter, apower, bexpn, bcenter, bpower )

    overlapx1 = overlap( aexpn, acenter, apower, bexpn, bcenter, (bI+2,bJ,bK) )

    overlapy1 = overlap( aexpn, acenter, apower, bexpn, bcenter, (bI,bJ+2,bK) )

    overlapz1 = overlap( aexpn, acenter, apower, bexpn, bcenter, (bI,bJ,bK+2) )

    overlapx2 = overlap( aexpn, acenter, apower, bexpn, bcenter, (bI-2,bJ,bK) )

    overlapy2 = overlap( aexpn, acenter, apower, bexpn, bcenter, (bI,bJ-2,bK) )

    overlapz2 = overlap( aexpn, acenter, apower, bexpn, bcenter, (bI,bJ,bK-2) )

    term0 = bexpn*( 2 * (bI + bJ + bK) + 3)*overlap0
    term1 = -2*(bexpn^2)*(overlapx1 + overlapy1 + overlapz1)
    term2 = -0.5*(bI*(bI-1)*overlapx2+bJ*(bJ-1)*overlapy2 + bK*(bK-1)*overlapz2)
    return term0 + term1 + term2
end

kinetic( a::CGBF, b::CGBF ) = contract( kinetic, a, b )


text:
using Documenter, Pumas
using Bioequivalence

ENV["COLUMNS"] = 120
ENV["LINES"] = 30

DocMeta.setdocmeta!(Pumas, :DocTestSetup, :(using Pumas, Bioequivalence); recursive=true)
DocMeta.setdocmeta!(Bioequivalence, :DocTestSetup, :(using Bioequivalence, Random, DataFrames); recursive=true)

makedocs(
  modules=Module[Pumas, Bioequivalence],
  doctest=true,
  clean=true,
  format =Documenter.HTML(assets = [joinpath("assets", "custom.css")]),
  sitename="Pumas",
  authors="PumasAI",
  pages = Any[
    "Home" => "index.md",
    "Tutorials" => Any[
      "tutorials/introduction.md",
    ],
    "Basics" => Any[
      "basics/overview.md",
      "basics/models.md",
      "basics/doses_subjects_populations.md",
      "basics/simulation.md",
      "basics/estimation.md",
      "basics/nca.md",
      "basics/be.md",
      "basics/faq.md",
    ],
    "Model Components" => Any[
      "model_components/domains.md",
      "model_components/dosing_control.md",
      "model_components/dynamical_types.md",
    ],
    "Diagnostics" => Any[
      "analysis/diagnostics.md",
    ],
  ]
  )

deploydocs(
   repo = "github.com/PumasAI/PumasDocs.jl.git",
   push_preview = true,
)


text:
#=
using Test
using LatticeTools

struct Phase{T<:Number} <: AbstractSymmetryOperation
    value::T
    function Phase(value::T) where {T<:Number}
        if !isapprox(abs(value), one(value))
            throw(ArgumentError("value $value not of unit length"))
        end
        return new{T}(value)
    end
end

LatticeTools.isidentity(arg::Phase) = isone(arg)
Base.:(*)(lhs::Phase, rhs::Phase) = Phase(lhs.value * rhs.value)
Base.:(^)(lhs::Phase, rhs::Real) = Phase(lhs.value^rhs)
Base.inv(lhs::Phase) = Phase(conj(lhs.value))
Base.isapprox(lhs::Phase, rhs::Phase; kwargs...) = isapprox(lhs.value, rhs.value; kwargs...)
Base.hash(p::P, h::UInt) where {P<:Phase} = Base.hash(P, Base.hash(p.value, h))

@testset "DirectProductOperation" begin
    α = Phase(-1)
    θ = Phase(cis(2π/3))
    t = TranslationOperation([1, 0])
    let
        θt = directproduct(θ, t)
        θt2 = directproduct(Phase(cis(2π/3)), TranslationOperation([1, 0]))
        # @show θt
        # @show θt2
        # @show (θt.operations) == (θt2.operations)
        # @show typeof(θt.operations)
        # @show typeof(θt2.operations)
        # @show hash(θt.operations), hash(θt2.operations)
        # @show (θt.operations), (θt2.operations)

        @test θt == θt2
        @test hash(θt) == hash(θt2)
        @test θt != t ×ˢ θ
        @test θt*θt == directproduct(θ*θ, t*t)
        @test θt^3 == directproduct(θ^3, t^3)
        @test inv(θt) == directproduct(inv(θ), inv(t))
    end
    let
        θt = θ ×ˢ t
        @test θt^3 == directproduct(θ^3, t^3)
        @test inv(θt) == directproduct(inv(θ), inv(t))
    end
    let
        αθt = (α ×ˢ θ) ×ˢ t
        @test αθt == α ×ˢ (θ ×ˢ t)
        @test αθt == directproduct(directproduct(α, θ), t)
        @test αθt == directproduct(α, directproduct(θ, t))
        @test αθt == directproduct(α, directproduct(θ, t))
    end
    @test directproduct([α, α^2], [t, t^2, t^3]) == [
        (α ×ˢ t)   (α ×ˢ t^2)   (α ×ˢ t^3) ;
        (α^2 ×ˢ t) (α^2 ×ˢ t^2) (α^2 ×ˢ t^3)
    ]
    @test [α, α^2] ×ˢ [t, t^2, t^3] == [
        (α ×ˢ t)   (α ×ˢ t^2)   (α ×ˢ t^3) ;
        (α^2 ×ˢ t) (α^2 ×ˢ t^2) (α^2 ×ˢ t^3)
    ]
end
=#


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Functions related to the state vector.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export satsv

"""
    function satsv(t::T1, r::AbstractVector{T2}, v::AbstractVector{T3} = [0,0,0], a::AbstractVector{T4} = [0,0,0]) where {T1<:Number, T2<:Number, T3<:Number, T4<:Number}
    function satsv(t::T1, vec::AbstractVector{T2}) where {T1<:Number, T2<:Number}

Create a new satellite state vector (see `SatelliteStateVector`) using the
position `r`, velocity `v`, and acceleration `a`. It is also possible to pass a
vector `vec` with the information concatenated.

!!! info

    The vectors `r`, `v`, and `a` must have at least 3 elements. In the case
    more elements are available, they will be neglected. On the other hand, the
    vector `v` must have 6 or 9 dimensions, indicating `[r;v]`, or `[r;v;a]`.

"""
function satsv(t::T1, r::AbstractVector{T2}, v::AbstractVector{T3},
               a::AbstractVector{T4} = [0,0,0]) where {T1<:Number, T2<:Number,
                                                       T3<:Number, T4<:Number}

    len_r = length(r)
    len_v = length(v)
    len_a = length(a)

    len_r  < 3 && error("The vector `r` must have at least 3 elements.")
    len_v  < 3 && error("The vector `v` must have at least 3 elements.")
    len_a  < 3 && error("The vector `a` must have at least 3 elements.")
    len_r != 3 && @warn("Only the first 3 elements of the vector `r` will be used!")
    len_v != 3 && @warn("Only the first 3 elements of the vector `v` will be used!")
    len_a != 3 && @warn("Only the first 3 elements of the vector `a` will be used!")

    T = promote_type(T1,T2,T3,T4)
    return SatelliteStateVector{T}(t = t, r = r[1:3], v = v[1:3], a = a[1:3])
end

function satsv(t::T1, vec::AbstractVector{T2}) where {T1<:Number, T2<:Number}
    len_vec = length(vec)

    len_vec ∉ [6,9] && error("The length of input vector must be 6, or 9.")

    if len_vec == 6
        return SatelliteStateVector(t = t, r = vec[1:3], v = vec[4:6])
    else
        return SatelliteStateVector(t = t, r = vec[1:3], v = vec[4:6], a = vec[7:9])
    end
end

################################################################################
#                                  Overloads
################################################################################

copy(sv::SatelliteStateVector) = SatelliteStateVector(sv.t, sv.r, sv.v, sv.a)

getindex(sv::SatelliteStateVector{T}, ::Colon) where T<:Number =
    SVector{9,T}(sv.r..., sv.v..., sv.a...)

function show(io::IO, sv::SatelliteStateVector{T}) where T
    t_str = @sprintf "%.g" sv.t

    print(io, "SatelliteStateVector{", string(T), "}: ")
    print(io, "t = ", t_str, " r = ", sv.r, " v = ", sv.v, " a = ", sv.a)
    return nothing
end

function show(io::IO, mime::MIME"text/plain", sv::SatelliteStateVector{T}) where T
    # Check if the `io` supports colors.
    color = get(io, :color, false)

    b = (color) ? _b : ""
    d = (color) ? _d : ""

    t_str = @sprintf "%.g" sv.t

    println(io, "SatelliteStateVector{", string(T), "}: ")
    println(io, "$b  t $d= ", t_str)
    println(io, "$b  r $d= ", sv.r)
    println(io, "$b  v $d= ", sv.v)
      print(io, "$b  a $d= ", sv.a)

    return nothing
end


text:
module SimpleCTF

export find_ctf
include("utils.jl")

using ArgParse
using LinearAlgebra: BLAS
using Base.Threads

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--pixelsize"
            help = "Pixel size in Angstroms"
            arg_type = Float64
            required = true
        "--voltage"
            help = "Accelerating voltage in kV"
            arg_type = Float64
            required = true
        "--Cs"
            help = "Spherical abberation in millimeters"
            arg_type = Float64
            required = true
        "--min_resolution"
            help = "Lower resolution for search"
            arg_type = Float64
            default = 25.
        "--max_resolution"
            help = "Upper resolution for search"
            arg_type = Float64
            default = 5.5
        "--amplitude_contrast"
            help = "Fraction of amplitude contrast"
            arg_type = Float64
            default = 0.07
        "--do_phase_search"
            help = "Additionally search for phase shift"
            action = :store_true
        "--csv"
            help = "csv to store results"
            default = "simplectf_results.csv"
        "images"
            nargs = '*'
            default = []
            help = "Aligned images for ctf estimation. Must be tif format."
            required = true
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()

    pixelsize = parsed_args["pixelsize"]
    voltage = parsed_args["voltage"]
    Cs = parsed_args["Cs"]
    min_resolution = parsed_args["min_resolution"]
    max_resolution = parsed_args["max_resolution"]
    amplitude_contrast = parsed_args["amplitude_contrast"]
    do_phase_search = parsed_args["do_phase_search"]
    out_csv = parsed_args["csv"]
    input_images = parsed_args["images"]

    BLAS.set_num_threads(1)
    FFTW.set_num_threads(1)

    results = []
    @time @threads for input_image in input_images
        println("processing $input_image")
        reference_spectrum, estimated_spectrum, search_result = find_ctf(
            input_image,
            pixelsize,
            voltage,
            Cs,
            min_resolution,
            max_resolution,
            amplitude_contrast,
            do_phase_search,
        )
        push!(results, (input_image, search_result))
        #diagnostic = [reference_spectrum[end:-1:1, :]; estimated_spectrum]
        #diagnostic ./= maximum(diagnostic)
    end

    f = open(out_csv, "w")
    write(f, "filename,defocus1,defocus2,astigmatism angle,phase shift,correlation\n")
    for (input_image, search_result) in results
        write(f, "$input_image,$search_result\n")
    end
    close(f)
    println("results written to $out_csv")

end

end


text:
type DeploymentProblem{ IM <: AbstractMatrix{Int},
                        BM <: AbstractMatrix{Bool}}
    nambulances::Int
    nlocations::Int
    nregions::Int
    demand::IM      # nperiod x nregion
    train::Vector{Int}
    test::Vector{Int}
    coverage::BM    # nregion x nlocation
    adjacency::BM   # nregion x nregion
end

function DeploymentProblem{BM}(
        hourly_calls::DataFrame,
        adjacent_nbhd::DataFrame,
        coverage::BM;
        namb = 30,
        train_filter = (hourly_calls[:year] .== 2012) .* (hourly_calls[:month] .<= 3)
    )
    regions = Int[parse(Int,string(x)[2:end]) for x in names(hourly_calls[5:end])]
    locations = collect(1:size(coverage,2))
    adjacent = convert(Array, adjacent_nbhd[2:end])[regions,regions] .> 0.5
    demand = convert(Array,hourly_calls[:,5:end])

    indices = 1:nrow(hourly_calls)
    train_indices = indices[train_filter]
    test_indices = indices[!train_filter]

    DeploymentProblem(
        namb,
        length(locations),
        length(regions),
        demand,
        train_indices,
        test_indices,
        coverage[regions,:],
        Array{Bool,2}(adjacent)
    )
end

function naive_solution(p::DeploymentProblem)
    # evenly distribute the ambulances over all the locations
    x = zeros(Int, p.nlocations)
    for i in 0:p.nambulances-1
        x[1 + (i % p.nlocations)] += 1
    end
    x
end

type DispatchProblem
    emergency_calls::DataFrame
    hospitals::DataFrame
    stations::DataFrame
    coverage::Matrix{Bool} # (nbhd x stns)
    turnaround::LogNormal
    deployment::Vector{Int}

    wait_queue::Vector{Vector{Int}} # length nbhd
    available::Vector{Int} # length stns
    redeploy_events::Vector{Tuple{Int,Int,Int,Int}} # amb,from,to,time

    DispatchProblem{BM}(emergency_data::DataFrame,
                        hospitals::DataFrame,
                        stations::DataFrame,
                        coverage::BM;
                        turnaround::LogNormal = LogNormal(3.65, 0.3)) =
        new(emergency_data, hospitals, stations, coverage, turnaround)
end

function initialize!(problem::DispatchProblem, deployment::Vector{Int})
    problem.wait_queue = [Int[] for i in 1:size(problem.coverage,1)]
    problem.available = copy(deployment)
    problem.deployment = deepcopy(deployment)
    problem.redeploy_events = Tuple{Int,Int,Int,Int}[]

    problem.emergency_calls[:arrival_seconds] =
        cumsum(problem.emergency_calls[:interarrival_seconds])

    problem
end

function DispatchProblem{BM}(
        emergency_data::DataFrame,
        hospitals::DataFrame,
        stations::DataFrame,
        coverage::BM,
        deployment::Vector{Int};
        turnaround::LogNormal = LogNormal(3.65, 0.3)
    )
    problem = DispatchProblem(
        emergency_data,
        hospitals,
        stations,
        coverage,
        turnaround=turnaround
    )
    initialize!(problem, deployment)
end

function returned_to!(problem::DispatchProblem, location::Int, t::Int)
    @assert problem.available[location] >= 0
    problem.available[location] += 1
end

type Params
    α::Float64 # Probabilistic Guarantee
    ε::Float64 # Convergence
    δ::Float64 # Solver Tolerance

    nperiods::Int # for StochasticDeployment

    maxiter::Int # for RobustDeployment
end
Params() = Params(0.01, 0.5, 1e-6, 500, 50)
params = Params()


text:
# This file defines the BaseNode type. This type defines a network node that a #   personnel member can be in.

export BaseNode
"""
The `BaseNode` type defines a node that a personnel member in the simulation can have. A node in this instance is a collection of attributes and values that must be satisfied before a personnel member is considered to be in that node.

The type contains the following fields:
* `name::String`: the name of the node.
* `target::Int`: the target number of personnel members in this node. Default = -1 (no target)
* `attrition::String`: The name of the attrition scheme attached to this node. Default = "default"
* `requirements::Dict{String,String}`: the attributes (key) and the value(s) these attributes must have (value) for the personnel member to have this particular node.

The type also contains one other field, used only during the simulation:
* `inNodeSince::Dict{String,Float64}`: a record of all the persons in the   node, including the time they last entered the node.

Constructor:
```
BaseNode( name::String )
```
This constructor generates a `BaseNode` with name `name`.
"""
mutable struct BaseNode

    name::String
    target::Int
    attrition::String
    requirements::Dict{String,String}

    inNodeSince::Dict{String,Float64}

    # Basic constructor.
    function BaseNode( name::String )::BaseNode

        newNode = new()
        newNode.name = name
        newNode.target = -1
        newNode.attrition = "default"
        newNode.requirements = Dict{String,String}()

        newNode.inNodeSince = Dict{String,Float64}()
        return newNode

    end  # BaseNode( name )

end  # type BaseNode


text:
"""
# Method p47

Analysis of plates (Plane structural element) using 4-node Quadrilateral
finite elements. Homogeneous material with identical elements. Mesh
numbered in x or y direction.

### Constructors
```julia
p47(data)
```
### Arguments
```julia
* `data::Dict{Symbol, Any}` : Dictionary containing all input data
```

### Required data dictionary keys
```julia
* struc_el::StructuralElement                          : Structural element
* support::Array{Tuple{Int,Array{Int,1}},1}        : Fixed-displacements vector
* loaded_nodes::Array{Tuple{Int,Array{Float64,1}},1} : Node load vector
* properties::Vector{Float64}                          : Material properties
* x_coords::FloatRange{Floalt64}                       : x-coordinate vector
* y_coords::FloatRange{Floalt64}                       : y-coordinate vector
* thickness:: Float64                                  : Thickness of plate
```

### Optional additional data dictionary keys
```julia
* penalty = 1e20               : Penalty used for fixed degrees of freedoms
* etype::Vector{Int}         : Element material vector if np_types > 1
```

### Return values
```julia
* (fm_df, sigma_df)            : Tuple of jFem, dis_df and fm_df
                                  where:
                                    fm_df         : Forces and moments data table
                                    sigma_df      : Stresses data table
```

### Related help
```julia
?StructuralElement             : List of available structural element types
?Plane                         : Help on a Plane structural element
?FiniteElement                 : List finite element types
?Quadrilateral                 : Help on Quadrilateral finite element
```
"""
function p47(data::Dict{Symbol, Any})

  # Setup basic dimensions of arrays

  # Parse & check FEdict data

  if :struc_el in keys(data)
    struc_el = data[:struc_el]
  else
    println("No fin_el type specified.")
    return
  end

  ndim = struc_el.ndim
  nst = struc_el.nst

  # Add radial stress
  if ndim == 3 && struc_el.axisymmetric
    nst = 4
  end

  fin_el = struc_el.fin_el
  @assert typeof(fin_el) <: FiniteElement

  if typeof(fin_el) == Line
    (nels, nn) = mesh_size(fin_el, struc_el.nxe)
  elseif typeof(fin_el) == Triangle || typeof(fin_el) == Quadrilateral
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye)
  elseif typeof(fin_el) == Hexahedron
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye, struc_el.nze)
  else
    println("$(typeof(fin_el)) is not a known finite element.")
    return
  end

  nodof = fin_el.nodof           # Degrees of freedom per node
  ndof = fin_el.nod * nodof      # Degrees of freedom per fin_el

  # Update penalty if specified in FEdict

  penalty = 1e20
  if :penalty in keys(data)
    penalty = data[:penalty]
  end

  # Allocate all arrays

  # Start with arrays to be initialized from FEdict

  if :properties in keys(data)
    prop = zeros(size(data[:properties], 1), size(data[:properties], 2))
    for i in 1:size(data[:properties], 1)
      prop[i, :] = data[:properties][i, :]
    end
  else
    println("No :properties key found in FEdict")
  end

  nf = ones(Int, nodof, nn)
  if :support in keys(data)
    for i in 1:size(data[:support], 1)
      nf[:, data[:support][i][1]] = data[:support][i][2]
    end
  end

  x_coords = zeros(nn)
  if :x_coords in keys(data)
    x_coords = data[:x_coords]
  end

  y_coords = zeros(nn)
  if :y_coords in keys(data)
    y_coords = data[:y_coords]
  end

  z_coords = zeros(nn)
  if :z_coords in keys(data)
    z_coords = data[:z_coords]
  end

  etype = ones(Int, nels)
  if :etype in keys(data)
    etype = data[:etype]
  end

  # All other arrays

  points = zeros(struc_el.nip, ndim)
  g = zeros(Int, ndof)
  g_coord = zeros(ndim,nn)
  fun = zeros(fin_el.nod)
  coord = zeros(fin_el.nod, ndim)
  gamma = zeros(nels)
  jac = zeros(ndim, ndim)
  g_num = zeros(Int, fin_el.nod, nels)
  der = zeros(ndim, fin_el.nod)
  deriv = zeros(ndim, fin_el.nod)
  bee = zeros(nst,ndof)
  km = zeros(ndof, ndof)
  mm = zeros(ndof, ndof)
  gm = zeros(ndof, ndof)
  kg = zeros(ndof, ndof)
  eld = zeros(ndof)
  weights = zeros(struc_el.nip)
  g_g = zeros(Int, ndof, nels)
  num = zeros(Int, fin_el.nod)
  actions = zeros(ndof, nels)
  displacements = zeros(size(nf, 1), ndim)
  gc = ones(ndim)
  dee = zeros(nst,nst)
  sigma = zeros(nst)
  axial = zeros(nels)

  formnf!(nodof, nn, nf)
  neq = maximum(nf)

  # Find global array sizes

  for iel in 1:nels
    geom_rect!(fin_el, iel, x_coords, y_coords, coord, num, struc_el.direction)
    num_to_g!(num, nf, g)
    g_num[:, iel] = num
    g_coord[:, num] = coord'
    g_g[:, iel] = g
  end
  println()

  println("There are $(neq) equations.")

  sample!(fin_el, points, weights)

  d2x = zeros(neq)
  d2y = zeros(neq)
  d2xy = zeros(neq)
  dtd = zeros(ndof, ndof)
  @assert :thickness in keys(data)
  aa = x_coords[2]-x_coords[1]
  bb = y_coords[2]-y_coords[1]

  gsm = spzeros(neq, neq)
  for iel in 1:nels
    e = prop[etype[iel], 1]
    v = prop[etype[iel], 2]
    d = e * data[:thickness]^3 / (12.0 * (1.0 - v * v))
    g = g_g[:, iel]
    fill!(km, 0)
    for i in 1:struc_el.nip
      fmplat!(d2x, d2y, d2xy, points, aa, bb, i)
      for k in 1:ndof
        tmp = (d2x[k]*d2x[:]/aa^4 + d2y[k]*d2y[:]/bb^4 +
        (v*d2x[k]*d2y[:]+v*d2x[:]*d2y[k] + 2.0*(1.0-v)*d2xy[k]*d2xy[:])/(aa^2*bb^2))
        dtd[k, :] = 4.0aa*bb*d*weights[i]*tmp
        dtd[:, k] = dtd[k, :]
      end
      km = km + dtd
    end
    fsparm!(gsm, g, km)
  end
  println()

  loads = OffsetArray(zeros(neq+1), 0:neq)
  if :loaded_nodes in keys(data)
    for i in 1:size(data[:loaded_nodes], 1)
      loads[nf[:, data[:loaded_nodes][i][1]]] = data[:loaded_nodes][i][2]
    end
  end

  fixed_freedoms = 0
  if :fixed_freedoms in keys(data)
    fixed_freedoms = size(data[:fixed_freedoms], 1)
  end
  no = zeros(Int, fixed_freedoms)
  node = zeros(Int, fixed_freedoms)
  sense = zeros(Int, fixed_freedoms)
  value = zeros(Float64, fixed_freedoms)
  if :fixed_freedoms in keys(data) && fixed_freedoms > 0
    for i in 1:fixed_freedoms
      no[i] = nf[data[:fixed_freedoms][i][2], data[:fixed_freedoms][i][1]]
      value[i] = data[:fixed_freedoms][i][3]
      gsm[no[i], no[i]] += penalty
      loads[no[i]] = gsm[no[i], no[i]] * value
    end
  end

  cfgsm = cholesky(gsm)
  loads[1:neq] = cfgsm \ loads[1:neq]

  disp = Float64[]
  rotx = Float64[]
  roty = Float64[]
  twistxy = Float64[]
  for i in 1:nn
    append!(disp, loads[nf[1,i]])
    append!(rotx, loads[nf[2,i]])
    append!(roty, loads[nf[3,i]])
    append!(twistxy, loads[nf[4,i]])
  end
  dis_df = DataFrame(
    disp = disp,
    rotx = rotx,
    roty = roty,
    twistxy = twistxy
  )

  struc_el.nip = 1
  points = zeros(struc_el.nip, ndim)
  weights = zeros(struc_el.nip)

  sample!(fin_el, points, weights)

  sigx = Float64[]
  sigy = Float64[]
  tauxy = Float64[]

  bm = Vector{Float64}
  for iel in 1:nels
    e = prop[etype[iel], 1]
    v = prop[etype[iel], 2]
    d = e * data[:thickness]^3 / (12.0 * (1.0 - v * v))
    g = g_g[:, iel]
    for i in 1:struc_el.nip
      fmplat!(d2x, d2y, d2xy, points, aa, bb, i)
      bm = zeros(3)
      for k in 1:ndof
        bm[1] += 4.0*d*(d2x[k]/aa/aa + v*d2y[k]/bb/bb)*loads[g[k]]
        bm[2] += 4.0*d*(v*d2x[k]/aa/aa + d2y[k]/bb/bb)*loads[g[k]]
        bm[3] += 4.0*d*(1.0-v)*(d2xy[k]/aa/bb*loads[g[k]])
      end
      append!(sigx, bm[1])
      append!(sigy, bm[2])
      append!(tauxy, bm[3])
    end
  end

  sigma_df = DataFrame(
    sigx = sigx,
    sigy = sigy,
    tauxy = tauxy
  )

  (dis_df, sigma_df)
end

max_stars_repo_path: src/4 Static Equilibrium/p47.jl

max_issues_repo_path: src/4 Static Equilibrium/p47.jl

max_forks_repo_path: src/4 Static Equilibrium/p47.jl


text:
#  Unit "Jade.jl" of the Diagonalization.jl Package for Julia language
#
#  MIT License
#  Copyright (c) 2019, 2020
#  Marco Congedo, CNRS, Grenoble, France:
#  https://sites.google.com/site/marcocongedo/home
#  Konstantin Usevich, CNRS, Nancy, France
#  http://w3.cran.univ-lorraine.fr/konstantin.usevich/?q=content/home

# ? CONTENTS :
#  This unit implements two orthogonal AJD algorithms:
#  1) The famous JADE algorithm of Cardoso and Souloumiac (1996),
#  performing GIVENS ROTATIONS of p,q pairs of vectors according to
#  the CYCLIC SCHEME, i.e., all p,q pairs are rotated in the order
#  for all p=1:n-1, q=i+1:n.
#  The code here below is optimized for both real and complex data input and
#  has been adapted in Julia from code freely made available from the authors.
#  J.-F. Cardoso, A. Souloumiac (1996) Jacobi angles for simultaneous
#  diagonalization. SIAM Journal on Matrix Analysis and Applications,
#  17(1), 161–164.
#  2) The JADEmax algorithm proposed by Usevich, Li and Comon (2020).
#  This uses the same GIVENS ROTATIONS however the p,q pairs are chosen
#  as the pair corresponding to the maximum p,q entry of the Riemannian
#  gradient of the cost function. The convergence is much faster in this case
#  at an extra cost of updating the gradient after each rotation.
#  Simulations and tests on real data indicate that JADEmax is overall faster
#  for n>20.
#
#  NB: JADE may be optimized by multi-threading the optimization of the pairs
#  e.g., using the round-Robin tournament scheme.
#  On the other hand, JADEmax is by definition a sequential procedure.

#  These algorithms handle the AJD diagonalization procedure corresponding
#  to the case m=1, k>1 according to the taxonomy adopted in this package.

## Compute the Givens angle (sine, cosine) for p,q pair
# and corrsponding p,q gradient element for real data
@inline function givensAngle(C::Matrix{T}, p, q, 𝓹, 𝓺,
							  e₁, e₂, e₃, Φ, Φₜ) where T <:Real
  e₁[:] = C[p, 𝓹] - C[q, 𝓺]
  e₂[:] = C[p, 𝓺] + C[q, 𝓹]
  a = e₁⋅e₁ - e₂⋅e₂
  b = 2. * e₁⋅e₂
  θ = 0.5 * atan(b, a + √(a^2 + b^2))
  s, c = sincos(θ)
  return c, s, abs(s), abs(b/4) # b/4 is the gradient
end

# Compute the Givens angle (sine, cosine) for p,q pair
# and corrsponding p,q gradient element for complex data
@inline function givensAngle(C::Matrix{T}, p, q, 𝓹, 𝓺,
								e₁, e₂, e₃, Φ, Φₜ) where T <:Complex
  e₁[:] = C[p, 𝓹] - C[q, 𝓺]
  e₂[:] = C[p, 𝓺]
  e₃[:] = C[q, 𝓹]
  E = Hermitian(UpperTriangular{T}([e₁⋅e₁ e₂⋅e₁ e₃⋅e₁; 0. e₂⋅e₂ e₃⋅e₂; 0. 0. e₃⋅e₃]))
  Γ = real(Φ*E*Φₜ)
  𝛉 = eigvecs(Γ)[:, 3] # Julia sorts the eigvecs; the 3rd<->max ev
  if real(𝛉[1])<0. 𝛉 = -𝛉 end  # added 'real' to prevent error
  c = √(0.5 + 𝛉[1]/2.)         # cosine
  s = 0.5*(𝛉[2] - 𝛉[3]im)/c   # sine
  return c, s, abs(s), abs(Γ[1, 2]+Γ[1, 3]im) # Γ[1, 2]+Γ[1, 3]im is the gradient
end

# Update U and all Cₖ matrices by a Givens rotation (both real and complex data)
@inline function jadeUpdate!(U, C, c, s, p, q, 𝓹, 𝓺)
	G = [c -conj(s); s c]
	⊶ = [p, q]  # p,q index pair
	U[:, ⊶] = U[:, ⊶]*G
	C[⊶, :] = G'*C[⊶, :]
	C[:, [𝓹 𝓺]] = [c*C[:, 𝓹]+s*C[:, 𝓺] -conj(s)*C[:, 𝓹]+c*C[:, 𝓺]]
end

# Do one cyclic sweep: rotate p,q pairs for all p=1:n-1, q=i+1:n
# For each p,q pair update U and all Cₖ if max abs(sine(p,q angle))>tol
# Return max abs(sine(p,q angle)) over all pairs
@inline function cyclicSweepAngle!(C, e₁, e₂, e₃, Φ, Φₜ, U, n, nk, tolerance)
  ∡ = 0.
  @inbounds for p = 1:n-1
	 𝓹 = p:n:nk
	 for q = p+1:n
		𝓺 = q:n:nk
		c, s, 𝓈, g = givensAngle(C, p, q, 𝓹, 𝓺, e₁, e₂, e₃, Φ, Φₜ)
		𝓈 > tolerance ? jadeUpdate!(U, C, c, s, p, q, 𝓹, 𝓺) : nothing
		∡ = max(∡, 𝓈)  # 𝓈 is abs(sine of the p,q Givens angle)
	 end
  end
  return ∡
end

# Do one cyclic sweep: rotate p,q pairs for all p=1:n-1, q=i+1:n
# For each p,q pair update U and all Cₖ if abs(p,q gradient's element)>tol
# Return the maximum of abs(p,q gradient's element) over all pairs
@inline function cyclicSweepGradient!(C, e₁, e₂, e₃, Φ, Φₜ, U, n, nk, tolerance)
  ∇pq = 0.
  @inbounds for p = 1:n-1
	 𝓹 = p:n:nk
	 for q = p+1:n
		𝓺 = q:n:nk
		c, s, 𝓈, g = givensAngle(C, p, q, 𝓹, 𝓺, e₁, e₂, e₃, Φ, Φₜ)
		g > tolerance ? jadeUpdate!(U, C, c, s, p, q, 𝓹, 𝓺) : nothing
		∇pq = max(∇pq, g)  # g is abs(of p,q Riemannian gradient's element)
	 end
  end
  return ∇pq
end

# Compute and return the Riemannian gradient (real data)
@inline function jadeRiemannian∇(C::Matrix{T}, D, n, k, 𝓹₁, 𝓹₂) where T <:Real
	∇=zeros(T, n, n)
	@inbounds @simd for i=1:k
		Cₖ = C[:, 𝓹₁[i]: 𝓹₂[i]]
		D[i][:] = diag(Cₖ)
		∇ += Cₖ.*D[i]'
	end
	return ∇-∇'
end

# Compute and return the Riemannian gradient (real data)
@inline function jadeRiemannian∇(C::Matrix{T}, D, n, k, 𝓹₁, 𝓹₂) where T <:Complex
	∇=zeros(T, n, n)
	@inbounds @simd for i=1:k
		Cₖ = C[:, 𝓹₁[i]: 𝓹₂[i]]
		D[i][:] = conj(diag(Cₖ))
		∇ += Cₖ.*(transpose(D[i])) - Cₖ.*D[i] # the @. macro transform all ensuing operations in breadcasting element-wise operations
	end
	return ∇-∇'
end

# Update the p, q rows and columns of the gradient
@inline function jadeRiemannian∇Update!(∇, C::Matrix{T}, D, k, p, q, 𝓹, 𝓺, 𝓹₁, 𝓹₂) where T <:Real
	⊶ = [p, q]  # p,q index pair
	∇[:, ⊶] .= T(0)
	@inbounds @simd for i=1:k
		dp, dq = C[p, 𝓹[i]], C[q, 𝓺[i]]
		D[i][p], D[i][q] = dp, dq
		∇[:, ⊶] += [C[:, 𝓹[i]].*(dp.-D[i]) C[:, 𝓺[i]].*(dq.-D[i])]
	end
	∇[⊶, :]=-∇[:, ⊶]'
end

# Update the p, q rows and columns of the gradient
@inline function jadeRiemannian∇Update!(∇, C::Matrix{T}, D, k, p, q, 𝓹, 𝓺, 𝓹₁, 𝓹₂) where T <:Complex
	⊶ = [p, q]  # p,q index pair
	∇[:, ⊶] .= T(0)
	#@inbounds @simd for i=1:k
	for i=1:k
		dp, dq = conj(C[p, 𝓹[i]]), conj(C[q, 𝓺[i]])
		D[i][p], D[i][q] = dp, dq
		∇[:, ⊶] += [(C[:, 𝓹[i]].*(dp.-D[i]).+conj.(C[p, 𝓹₁[i]:𝓹₂[i]]).*conj.(dp.-D[i])) (C[:, 𝓺[i]].*(dq.-D[i]).+conj.(C[q, 𝓹₁[i]:𝓹₂[i]]).*conj.(dq.-D[i]))]
		#∇[:, ⊶] += [C[:, 𝓹[i]].*(dp-D[i])+conj(C[p, 𝓹₁[i]:𝓹₂[i]])*conj(dp.-D[i]) C[:, 𝓺[i]].*(dq.-D[i])+conj(C[q, 𝓹₁[i]:𝓹₂[i]])*conj(dq.-D[i])]
	end
	∇[⊶, :]=-∇[:, ⊶]'
end

# pre-allocate memory
@inline function jadeAllocateMemory(T, k)
	e₁ = zeros(T, k)
	e₂ = zeros(T, k)
	if T <:Complex
		e₃ = zeros(T, k)
		Φ = convert(Matrix{T}, [1 0 0; 0 1 1; 0 -im im])
		Φₜ = Φ'
	else
		e₃, Φ, Φₜ = nothing, nothing, nothing, nothing
	end
	return e₁, e₂, e₃, Φ, Φₜ
end

## PRIMITIVE JADE (Cyclic) algorithm:
#  It takes as input a n·nk matrix holding k horizontally stacked n·n real or
#  complex matrices, such as C=[C_1...C_k].
#  It finds an orthogoanl/unitary matrix U such that the
#  congruences U'*C_κ*U are as diagonal as possible for all κ=1:k.
#  `tol` is the convergence to be attained.
#  `maxiter` is the maximum number of iterations allowed.
#  if `updateRule`=`:angle` (default), the max abs of the sin of the angle
#  across all pairs is used as stopping criterion, otherwise (any other symbol)
#  the mas abs of the Riemannian gradient entry is used.
#  The former is the usual choice of the JADE algorithm, the latter allows to
#  compare directly the execution speed of JADE and JADEmax as JADEmax uses
#  that stopping criterion. Note that the Riemannian gradient criterion for J
#  JADE can be obtained at virtually no extra-cost from the computation of
#  the angle.
#  if `verbose`=true, the convergence attained at each iteration and other
#  information will be printed.
#
#  RETURN: B, the number of iterations and the convergence attained (a 3-tuple)
function jade(C::Matrix{T};
				tol 	 = 0.,
				maxiter = 1000,
				updateRule = :angle,
				verbose = false) where T<:Union{Real, Complex}

	(n, nk) = size(C)
	iter, conv, 😋, k, U = 1, 1., false, nk÷n, Matrix{T}(I, n, n)
	tolerance = tol≤0. ? √eps(real(T)) : tol
	e₁, e₂, e₃, Φ, Φₜ = jadeAllocateMemory(T, k)

	verbose && @info("Iterating JADE algorithm...")
	while true
		conv = updateRule==:angle ? cyclicSweepAngle!(C, e₁, e₂, e₃, Φ, Φₜ, U, n, nk, tolerance) :
			   						cyclicSweepGradient!(C, e₁, e₂, e₃, Φ, Φₜ, U, n, nk, tolerance)
		verbose && println("iteration: ", iter, "; convergence: ", conv)
	   	(overRun = iter == maxiter) && @warn("JADE reached the max number of iterations before convergence:", iter)
	   	(😋 = conv <= tolerance) || overRun==true ? break : iter += 1
	end
	verbose && @info("Convergence has "*(😋 ? "" : "not ")*"been attained.\n\n")

	return U, iter, conv
end

## PRIMITIVE JADEmax algorithm:
#  It takes as input a n·nk matrix holding k horizontally stacked n·n real or
#  complex matrices, such as C=[C_1...C_k].
#  It finds an orthogoanl/unitary matrix U such that the
#  congruences U'*C_κ*U are as diagonal as possible for all κ=1:k.
#  `tol` is the convergence to be attained.
#  `maxiter` is the maximum number of iterations allowed.
#  if `verbose`=true, the convergence attained at each iteration and other
#  information will be printed.
#
#  RETURN: B, the number of iterations and the convergence attained (a 3-tuple)
function jademax(C::Matrix{T};
					tol 	 = 0.,
					maxiter = 1000,
					verbose = false) where T<:Union{Real, Complex}
	(n, nk) = size(C)
	k, N, U, p, q = nk÷n, (n*(n-1))÷2, Matrix{T}(I, n, n), 0, 0
	𝓹₁, 𝓹₂ = 1:n:nk, n:n:nk
	𝓹, 𝓺 = 𝓹₁, 𝓹₂
	tolerance = tol≤0. ? √eps(real(T)) : tol
	iter, itN, itN0, conv, 😋 = 0, 1, 0, 1., false
	e₁, e₂, e₃, Φ, Φₜ = jadeAllocateMemory(T, k)
	D=[Vector{T}(undef, n) for i=1:k]

	verbose && @info("Iterating JADEmax algorithm...")
	# compute the full gradient only once at the beginning
	∇ = jadeRiemannian∇(C, D, n, k, 𝓹₁, 𝓹₂)
    while true
		# find max absolute element of the gradient and its position
		conv, maxpos = T <:Real ? findmax(∇) : findmax(abs2.(∇))
		p, q = maxpos[1], maxpos[2]
		𝓹, 𝓺  = p:n:nk, q:n:nk
		# find angle, update U, all Cₖ matrices and the Riemannian gradient
	  	c, s, ∡ = givensAngle(C, p, q, 𝓹, 𝓺, e₁, e₂, e₃, Φ, Φₜ)
		if conv > tolerance
			jadeUpdate!(U, C, c, s, p, q, 𝓹, 𝓺)
			jadeRiemannian∇Update!(∇, C, D, k, p, q, 𝓹, 𝓺, 𝓹₁, 𝓹₂)
		end

		(itN0=itN%N==0) && (iter+=1)
		verbose && itN0 && println("iteration(m): ", iter, "; convergence: ", conv)
		(overRun = iter == maxiter) && @warn("jademax reached the max number of iterations before convergence:", iter)
       	(😋 = conv <= tolerance) || overRun==true ? break : itN += 1
    end
	verbose && !itN0 && println("iteration(m): ", iter, " + ", itN%N, " sweeps ; convergence: ", conv)
    verbose && @info("Convergence has "*(😋 ? "" : "not ")*"been attained.\n\n")

	return U, iter, conv
end

#  ADVANCED JADE and JADEmax algorithms:
#  They take as input a vector of k real symmetric or complex Hermitian
#  matrices 𝐂 and finds an orthogoanl/Unitary matrix U such that the
#  congruences U'*𝐂_κ*U are as diagonal as possible for all κ=1:k.
#
#  if `trace1` is true (false by default), all input matrices are normalized
#  so as to have unit trace. Note that the diagonal elements
#  of the input matrices must be all positive.
#
#  `w` is an optional vector of k positive weights for each matrix in 𝐂.
#  if `w` is different from `nothing` (default), the input matrices are
#  weighted with these weights (after trace normalization if `trace1` is true).
#  A function can be passed as the `w` argument, in which case the kth weight
#  is found as the output of the function applied to the kth matrix in 𝐂.
#  A good choice in general is the `nonD` function declared in tools.jl unit.
#
#  if `whitening` = true is passed, the Arithmetic mean of the matrices in 𝐂 is
#  computed (using the PosDefManifold.jl package) and the matrices in 𝐂
#  are pre-transformed using the whitening matrix of the mean.
#  Dimensionality reduction can be obtained at this stage using optional
#  arguments `eVar` and `eVarMeth` (see documentation of the AJD constructors).
#
#  if sort=true (default) the column vectors of the B matrix are reordered
#  so as to sort in descending order the mean of the diagonal elements
#  of B'*𝐂_κ* over k.
#
#  if  `whitening` = false (default), a matrix can be provided with the `init`
#  argument in order to initialize B. In this case the actual AJD
#  will be given by init*B, where B is the output of the algorithms.
#
#  `tol` is the convergence to be attained.
#
#  `maxiter` is the maximum number of iterations allowed.
#
#  if `verbose`=true, the convergence attained at each iteration and other
#  information will be printed.
#
#  return: B, its pseudo-inverse, the mean diagonal elements of B'*mean(𝐂)*B,
#          the number of iterations and the convergence attained
function jade( 𝐂::Union{Vector{Hermitian}, Vector{Symmetric}};
               trace1   :: Bool  = false,
               w        :: Twf   = ○,
               preWhite :: Bool  = false,
               sort     :: Bool  = true,
               init     :: Union{Matrix, Nothing} = ○,
               tol      :: Real  = 0.,
               maxiter  :: Int   = 1000,
               verbose  :: Bool  = false,
            eVar     :: TeVaro = ○,
            eVarMeth :: Function = searchsortedfirst)

	# trace normalization and weighting
	𝐆 = _normalizeAndWeight(trace1, w, 𝐂)

	# pre-whiten or initialize and stack matrices horizontally
	W, C = _preWhiteOrInit(𝐂, preWhite, Euclidean, eVar, eVarMeth, init, :stacked)

	(n, nk) = size(C)

	U, iter, conv = jade(C; tol=tol, maxiter=maxiter, verbose=verbose)

	# permute the vectors of U
	D=Diagonal([mean(C[i, i:n:nk]) for i=1:n])
	λ = sort ? _permute!(U, D, n) : diag(D)

	return preWhite ? (W.F*U, U'*W.iF, λ, iter, conv) :
	                  (U, Matrix(U'), λ, iter, conv)
end

function jademax( 𝐂::Union{Vector{Hermitian}, Vector{Symmetric}};
               trace1   :: Bool  = false,
               w        :: Twf   = ○,
               preWhite :: Bool  = false,
               sort     :: Bool  = true,
               init     :: Union{Matrix, Nothing} = ○,
               tol      :: Real  = 0.,
               maxiter  :: Int   = 1000,
               verbose  :: Bool  = false,
            eVar     :: TeVaro = ○,
            eVarMeth :: Function = searchsortedfirst)

	# trace normalization and weighting
	𝐆 = _normalizeAndWeight(trace1, w, 𝐂)

	# pre-whiten or initialize and stack matrices horizontally
	W, C = _preWhiteOrInit(𝐂, preWhite, Euclidean, eVar, eVarMeth, init, :stacked)

	(n, nk) = size(C)

	U, iter, conv = jademax(C; tol=tol, maxiter=maxiter, verbose=verbose)

	# permute the vectors of U
	D=Diagonal([mean(C[i, i:n:nk]) for i=1:n])
	λ = sort ? _permute!(U, D, n) : diag(D)

	return preWhite ? (W.F*U, U'*W.iF, λ, iter, conv) :
	                  (U, Matrix(U'), λ, iter, conv)
end


text:
import JuliaCLI

@async JuliaCLI.serve()  # TODO: don't
let path = joinpath(DEPOT_PATH[1], "jlcli", "socket")
    for i in 1:100
        isfile(path) && break
        sleep(0.1)
    end
end

include("cli_test.jl")


text:
export OrcJIT, OrcModule, OrcTargetAddress
export JITEventListener, GDBRegistrationListener, IntelJITEventListener,
       OProfileJITEventListener, PerfJITEventListener
export dispose, errormsg, compile!, remove!, add!,
       mangle, address, addressin, create_stub!, set_stub!,
       register!, unregister!, callback!
export JITTargetMachine

@checked struct OrcJIT
    ref::API.LLVMOrcJITStackRef
end

Base.unsafe_convert(::Type{API.LLVMOrcJITStackRef}, orc::OrcJIT) = orc.ref

"""
    OrcJIT(::TargetMachine)

Creates a OrcJIT stack based on the provided target machine.

!!! warning
    Takes ownership of the provided target machine.
"""
function OrcJIT(tm::TargetMachine)
    OrcJIT(API.LLVMOrcCreateInstance(tm))
end

function dispose(orc::OrcJIT)
    API.LLVMOrcDisposeInstance(orc)
end

function errormsg(orc::OrcJIT)
    # The error message is owned by `orc`, and will
    # be disposed along-side the OrcJIT.
    unsafe_string(LLVM.API.LLVMOrcGetErrorMsg(orc))
end

struct OrcModule
    handle::API.LLVMOrcModuleHandle
end
Base.convert(::Type{API.LLVMOrcModuleHandle}, mod::OrcModule) = mod.handle

function compile!(orc::OrcJIT, mod::Module, resolver = C_NULL, ctx = C_NULL; lazy=false)
    r_mod = Ref{API.LLVMOrcModuleHandle}()
    if lazy
        API.LLVMOrcAddLazilyCompiledIR(orc, r_mod, mod, resolver, ctx)
    else
        API.LLVMOrcAddEagerlyCompiledIR(orc, r_mod, mod, resolver, ctx)
    end
    OrcModule(r_mod[])
end

function Base.delete!(orc::OrcJIT, mod::OrcModule)
    LLVM.API.LLVMOrcRemoveModule(orc, mod)
end

function add!(orc::OrcJIT, obj::MemoryBuffer, resolver = C_NULL, ctx = C_NULL)
    r_mod = Ref{API.LLVMOrcModuleHandle}()
    API.LLVMOrcAddObjectFile(orc, r_mod, obj, resolver, ctx)
    return OrcModule(r_mod[])
end

function mangle(orc::OrcJIT, name)
    r_symbol = Ref{Cstring}()
    LLVM.API.LLVMOrcGetMangledSymbol(orc, r_symbol, name)
    symbol = unsafe_string(r_symbol[])
    LLVM.API.LLVMOrcDisposeMangledSymbol(r_symbol[])
    return symbol
end

struct OrcTargetAddress
    ptr::API.LLVMOrcTargetAddress
end
Base.convert(::Type{API.LLVMOrcTargetAddress}, addr::OrcTargetAddress) = addr.ptr

Base.pointer(addr::OrcTargetAddress) = reinterpret(Ptr{Cvoid}, addr.ptr % UInt) # LLVMOrcTargetAddress is UInt64 even on 32-bit

OrcTargetAddress(ptr::Ptr{Cvoid}) = OrcTargetAddress(reinterpret(UInt, ptr))

function create_stub!(orc::OrcJIT, name, initial)
    LLVM.API.LLVMOrcCreateIndirectStub(orc, name, initial)
end

function set_stub!(orc::OrcJIT, name, new)
    LLVM.API.LLVMOrcSetIndirectStubPointer(orc, name, new)
end

function address(orc::OrcJIT, name)
    r_address = Ref{API.LLVMOrcTargetAddress}()
    API.LLVMOrcGetSymbolAddress(orc, r_address, name)
    OrcTargetAddress(r_address[])
end

function addressin(orc::OrcJIT, mod::OrcModule, name)
    r_address = Ref{API.LLVMOrcTargetAddress}()
    API.LLVMOrcGetSymbolAddressIn(orc, r_address, mod, name)
    OrcTargetAddress(r_address[])
end

function callback!(orc::OrcJIT, callback, ctx)
    r_address = Ref{API.LLVMOrcTargetAddress}()
    LLVM.API.LLVMOrcLazyCompileCallback(orc, r_address, callback, ctx)
    return OrcTargetAddress(r_address[])
end

@checked struct JITEventListener
    ref::API.LLVMJITEventListenerRef
end
Base.unsafe_convert(::Type{API.LLVMJITEventListenerRef}, listener::JITEventListener) = listener.ref

function register!(orc::OrcJIT, listener::JITEventListener)
    LLVM.API.LLVMOrcRegisterJITEventListener(orc, listener)
end

function unregister!(orc::OrcJIT, listener::JITEventListener)
    LLVM.API.LLVMOrcUnregisterJITEventListener(orc, listener)
end

GDBRegistrationListener()  = JITEventListener(LLVM.API.LLVMCreateGDBRegistrationListener())
IntelJITEventListener()    = JITEventListener(LLVM.API.LLVMCreateIntelJITEventListener())
OProfileJITEventListener() = JITEventListener(LLVM.API.LLVMCreateOProfileJITEventListener())
PerfJITEventListener()     = JITEventListener(LLVM.API.LLVMCreatePerfJITEventListener())

function JITTargetMachine(;triple = LLVM.triple(),
                           cpu = "", features = "",
                           optlevel = LLVM.API.LLVMCodeGenLevelDefault)

    # Force ELF on windows,
    # Note: Without this call to normalize Orc get's confused
    #       and chooses the x86_64 SysV ABI on Win x64
    triple = LLVM.normalize(triple)
    if Sys.iswindows()
        triple *= "-elf"
    end
    target = LLVM.Target(triple=triple)
    @debug "Configuring OrcJIT with" triple cpu features optlevel

    tm = TargetMachine(target, triple, cpu, features,
                       optlevel,
                       LLVM.API.LLVMRelocStatic, # Generate simpler code for JIT
                       LLVM.API.LLVMCodeModelJITDefault, # Required to init TM as JIT
                       )
    return tm
end


text:
"""
Functionality for handling the default content of the various files (migrations, models, controllers, etc).
"""
module FileTemplates

import Inflector
using SearchLight

"""
    new_database_migration(module_name::String) :: String

Default content for a new SearchLight migration.
"""
function new_table_migration(module_name::String, resource::String) :: String
  Inflector.is_plural(resource) || (resource = Inflector.to_plural(resource))

  """
  module $module_name

  import SearchLight.Migrations: create_table, column, columns, pk, add_index, drop_table, add_indices

  function up()
    create_table(:$resource) do
      [
        pk()
        column(:column_name, :column_type)
        columns([
          :column_name => :column_type
        ])
      ]
    end

    add_index(:$resource, :column_name)
    add_indices($resource, :column_name_1, :column_name_2)
  end

  function down()
    drop_table(:$resource)
  end

  end
  """
end

"""
    new_relationship_table_migration(module_name::String) :: String

Default content for a new SearchLight migration.
"""
function new_relationship_table_migration(module_name::String, table_name::String, r1::String, r2::String) :: String
  """
  module $module_name

  import SearchLight.Migrations: create_table, column, columns, pk, add_index, drop_table

  function up()
    create_table(:$table_name) do
      [
        pk()
        column(:$(r1)_id, :int)
        column(:$(r2)_id, :int)
      ]
    end

    add_index(:$table_name, :$(r1)_id)
    add_index(:$table_name, :$(r2)_id)
  end

  function down()
    drop_table(:$table_name)
  end

  end
  """
end

function newmigration(module_name::String) :: String
  """
  module $module_name

  function up()

  end

  function down()

  end

  end
  """
end

"""
    newmodel(model_name::String, resource_name::String = model_name) :: String

Default content for a new SearchLight model.
"""
function newmodel(model_name::String, resource_name::String = model_name; pluralize::Bool = true) :: String
  """
  module $(Inflector.to_plural(model_name) |> uppercasefirst)

  import SearchLight: AbstractModel, DbId
  import Base: @kwdef

  export $model_name

  @kwdef mutable struct $model_name <: AbstractModel
    id::DbId = DbId()
  end

  end
  """
end

"""
    newvalidator(validator_name::String) :: String

Default content for a new SearchLight validator.
"""
function newvalidator(validator_name::String; pluralize::Bool = true) :: String
  """
  module $(Inflector.to_plural(validator_name) |> uppercasefirst)Validator

  using SearchLight, SearchLight.Validation

  function not_empty(field::Symbol, m::T, args::Vararg{Any})::ValidationResult where {T<:AbstractModel}
    isempty(getfield(m, field)) && return ValidationResult(invalid, :not_empty, "should not be empty")

    ValidationResult(valid)
  end

  function is_int(field::Symbol, m::T, args::Vararg{Any})::ValidationResult where {T<:AbstractModel}
    isa(getfield(m, field), Int) || return ValidationResult(invalid, :is_int, "should be an int")

    ValidationResult(valid)
  end

  function is_unique(field::Symbol, m::T, args::Vararg{Any})::ValidationResult where {T<:AbstractModel}
    findone(typeof(m); NamedTuple(field => getfield(m, field))... ) === nothing ||
      return ValidationResult(invalid, :is_unique, "already exists")

    ValidationResult(valid)
  end

  end
  """
end

function adapter_default_config end

function newconfig() :: String
  """
  env: $(SearchLight.config.app_env)

  $(adapter_default_config())
    config:
      log_queries: true
      log_level: :debug
  """
end

"""
    newtest(resource_name::String) :: String

Default content for a new test file.
"""
function newtest(resource_name::String) :: String
  """
  using Test, SearchLight, $(Inflector.to_plural(resource_name) |> uppercasefirst)

  @testset "$resource_name unit tests" begin

    ### Your tests here
    @test 1 == 1

  end;
  """
end

end


text:
function chebyshev_points(::Type{T}, s::Integer, ::Val{1}) where {T}
    shift_nodes(T[ @big sin( π*(s-2i+1) / (2s) ) for i in s:-1:1 ])
end

function chebyshev_points(::Type{T}, s::Integer, ::Val{2}) where {T}
    shift_nodes(T[ @big cos( π * (i-1) / (s-1) ) for i in s:-1:1 ])
end

"""
Gauss-Chebyshev quadrature.
"""
function GaussChebyshevQuadrature(::Type{T}, s::Integer) where {T}
    local tj::BigFloat = 0
    local th::BigFloat = 0

    c = chebyshev_points(BigFloat, s, Val(1))
    b = zero(c)
    for i in eachindex(b)
        tj = 0
        th = @big π * (2i-1) / (2s)
        for j in 1:div(s,2)
            tj += @big cos(2j*th) / (4j^2 - 1)
        end
        b[i] = @big (1 - 2tj) / s
    end

    QuadratureRule(2s, c, b, T)
end

GaussChebyshevQuadrature(s) = GaussChebyshevQuadrature(Float64, s)

"""
Lobatto-Chebyshev quadrature.
"""
function LobattoChebyshevQuadrature(::Type{T}, s::Integer) where {T}
    if s == 1
        throw(ErrorException("Lobatto-Chebyshev quadrature is not defined for one stage."))
    end

    local tj::BigFloat = 0
    local th::BigFloat = 0

    c = chebyshev_points(BigFloat, s, Val(2))
    b = zero(c)
    for i in eachindex(b)
        tj = 0
        th = @big π * i / (s+1)
        for j in 1:div(s+1,2)
            tj += @big sin((2j-1) * th) / (2j - 1)
        end
        b[i] = @big (2tj * sin(th)) / (s+1)
    end

    QuadratureRule(2s, c, b, T)
end

LobattoChebyshevQuadrature(s) = LobattoChebyshevQuadrature(Float64, s)

ChebyshevQuadrature(::Type{T}, s::Integer, ::Val{1}) where {T} = GaussChebyshevQuadrature(T,s)
ChebyshevQuadrature(::Type{T}, s::Integer, ::Val{2}) where {T} = LobattoChebyshevQuadrature(T,s)

ChebyshevQuadrature(s, kind) = ChebyshevQuadrature(Float64, s, Val(kind))


text:
import ACEbase: Discrete1pBasis
import ACE: evaluate, evaluate!

export Categorical1pBasis

# -------------------------

struct SList{N, T}
   list::SVector{N, T}

   function SList{N, T}(list::SVector{N, T})  where {N, T}
      if isabstracttype(T)
         error("`SList` can only contain a single type")
      end
      return new(list)
   end
end

SList(list::AbstractArray) = SList(SVector(list...))
SList(args...) = SList(SVector(args...))
SList(list::SVector{N, T}) where {N, T} = SList{N, T}(list)

Base.length(list::SList) = length(list.list)
Base.rand(list::SList) = rand(list.list)

i2val(list::SList, i::Integer) = list.list[i]

Base.iterate(list::SList, args...) = iterate(list.list, args...)

function val2i(list::SList, val)
   for j = 1:length(list)
      if list.list[j] == val
         return j
      end
   end
   error("val = $val not found in this list")
end

write_dict(list::SList{N,T}) where {N, T} =
      Dict( "__id__" => "ACE_SList",
                 "T" => write_dict(T),
              "list" => list.list )

function read_dict(::Val{:ACE_SList}, D::Dict)
   list = D["list"]
   T = read_dict(D["T"])
   svector = SVector{length(list), T}((T.(list))...)
   return SList(svector)
end

# -------------------------

@doc raw"""
`Categorical1pBasis` : defines the discrete 1p basis
```math
   \phi_q(u) = \delta(u - U_q),
```
where ``U_q, q = 1, \dots, Q`` are finitely many possible values that the
variable ``u`` may take. Suppose, e.g., we allow the values `[:a, :b, :c]`,
then
```julia
P = Categorical1pBasis([:a, :b, :c]; varsym = :u, idxsym = :q)
evaluate(P, State(u = :a))     # Bool[1, 0, 0]
evaluate(P, State(u = :b))     # Bool[0, 1, 0]
evaluate(P, State(u = :c))     # Bool[0, 0, 1]
```
If we evaluate it with an unknown state we get an error:
```julia
evaluate(P, State(u = :x))
# Error : val = x not found in this list
```
"""
struct Categorical1pBasis{VSYM, ISYM, LEN, T} <: Discrete1pBasis{Bool}
   categories::SList{LEN, T}
   label::String
end

_varsym(::Categorical1pBasis{VSYM, ISYM}) where {VSYM, ISYM} = VSYM
_isym(::Categorical1pBasis{VSYM, ISYM}) where {VSYM, ISYM} = ISYM

_val(X, B::Categorical1pBasis) = getproperty(X, _varsym(B))
_idx(b, B::Categorical1pBasis) = getproperty(b, _isym(B))

Base.length(B::Categorical1pBasis) = length(B.categories)

Categorical1pBasis(categories::AbstractArray;
              varsym::Symbol = nothing, idxsym::Symbol = nothing,
              label = "C$(idxsym)") =
      Categorical1pBasis(categories, varsym, idxsym, label)

Categorical1pBasis(categories::AbstractArray, varsym::Symbol, isym::Symbol, label::String) =
      Categorical1pBasis(SList(categories), varsym, isym, label)

Categorical1pBasis(categories::SList{LEN, T}, varsym::Symbol, isym::Symbol, label::String) where {LEN, T} =
      Categorical1pBasis{varsym, isym, LEN, T}(categories, label)

function evaluate(basis::Categorical1pBasis, X::AbstractState)
   A = Vector{Bool}(undef, length(basis))
   return evaluate!(A, basis, X)
end

function ACE.evaluate!(A, basis::Categorical1pBasis, X::AbstractState)
   fill!(A, false)
   A[val2i(basis.categories, _val(X, basis))] = true
   return A
end

symbols(basis::Categorical1pBasis) = [ _isym(basis), ]

indexrange(basis::Categorical1pBasis) = Dict( _isym(basis) => basis.categories.list )

isadmissible(b, basis::Categorical1pBasis) = (_idx(b, basis) in basis.categories)

get_index(B::Categorical1pBasis, b) = val2i(B.categories, _idx(b, B) )

degree(b, basis::Categorical1pBasis, args...) = 0

Base.rand(basis::Categorical1pBasis) = rand(basis.list)

function get_spec(basis::Categorical1pBasis, i)
   return NamedTuple{(_isym(basis),)}((i2val(basis.categories, i),))
end

get_spec(basis::Categorical1pBasis) = [ get_spec(basis, i) for i = 1:length(basis) ]

write_dict(B::Categorical1pBasis) =
      Dict( "__id__" => "ACE_Categorical1pBasis",
            "categories" => write_dict(B.categories),
            "VSYM" => String(_varsym(B)),
            "ISYM" => String(_isym(B)),
            "label" => B.label)

read_dict(::Val{:ACE_Categorical1pBasis}, D::Dict)  =
   Categorical1pBasis( read_dict(D["categories"]),
                  Symbol(D["VSYM"]), Symbol(D["ISYM"]),
                  D["label"] )


text:
import NeuroAnalysis
using Base.Test


text:
@testset "EpitopeMap(::String, ::Vector{String}, ::Vector{Int}, ::Vector{Int})" begin
    epitopes = ["SLYNTVATLY", "LYNTVATL"]
    start = [77, 78]
    stop = [86, 85]
    @test EpitopeMap("test", epitopes, start, stop,
        [parse_allele("B51"), parse_allele("A01")], 100) isa EpitopeMap

    @test_throws ErrorException EpitopeMap("", ["SLYNTVATL", ""], [77], [84],
        [parse_allele("A11")], 100) isa EpitopeMap
    @test_throws ErrorException EpitopeMap("", ["SLYNTVATLY"], [85], [77],
        [parse_allele("A03")], 100) isa EpitopeMap
    @test_throws ErrorException EpitopeMap("", ["SLYNTVATLY"], [-3], [5],
        [parse_allele("A01")], 100) isa EpitopeMap
    @test_throws ErrorException EpitopeMap("", ["SLYNTVATLY"], [77], [85],
        [parse_allele("A68")], 100) isa EpitopeMap
    @test_throws ErrorException EpitopeMap("test", epitopes, start, stop,
        [parse_allele("B51")], 40) isa EpitopeMap
end

@testset "epitope_map(::AbstractHLAData)" begin
    hla_types = [HLAType(parse_allele("A01", "A01", "B01", "B01", "C01", "C01")),
                 HLAType(parse_allele("A01", "A01", "B01", "B01", "C01", "C01")),
                 HLAType(parse_allele("A01", "A01", "B01", "B01", "C01", "C01")),
                 HLAType(parse_allele("A01", "A01", "B01", "B01", "C01", "C01")),
                 HLAType(parse_allele("A01", "A01", "B01", "B01", "C01", "C01"))]
    fasta_path = joinpath(@__DIR__, "data", "test.fasta")
    hla_data = HLAData("test", fasta_path, hla_types, missing, missing)

    map = @suppress Escape.epitope_map(hla_data)

    @test map isa Escape.EpitopeMap
    @test map.maplength == 14
end

@testset "Base.in(::Pair{Replacement, HLAAllele}, ::EpitopeMap)" begin
    epitopes = ["SLYNTVATLY", "LYNTVATL"]
    start = [77, 78]
    stop = [86, 85]
    alleles = [parse_allele("B51"), parse_allele("A01")]
    map = EpitopeMap("test", epitopes, start, stop, alleles, 100)
    r_1 = Replacement("test", 77, 'G', false)
    r_2 = Replacement("test", 76, 'S', false)
    r_3 = Replacement("test", 99, 'S', false)
    r_4 = Replacement("test", 80, 'G', false)

    @test (r_1 => parse_allele("B51")) in map
    @test !((r_1 => parse_allele("A01")) in map)
    @test !((r_1 => parse_allele("B52")) in map)
    @test !((r_2 => parse_allele("B51")) in map)
    @test !((r_2 => parse_allele("B51")) in map)
    @test !((r_3 => parse_allele("B51")) in map)
    @test !((r_3 => parse_allele("B52")) in map)
    @test ((r_4 => parse_allele("B51")) in map)
    @test ((r_4 => parse_allele("A01")) in map)
    @test !((r_4 => parse_allele("B52")) in map)

end

@testset "match_epitope(::String, ::HLAData)" begin
    data = Escape.HLADataset("Test").data[1]

    @test Escape.match_epitope("SARA", data) == 2
    @test Escape.match_epitope("SVMG", data) == 6
    @test isnothing(Escape.match_epitope("SMRA", data))
    @test Escape.match_epitope("SMRA", data, max_levenshtein_distance = 1) == 2
end


text:
using Documenter

if ispath("buld/")
    rm("build/", recursive=true)
end

if ispath("site/")
    rm("site/", recursive=true)
end

include("runWeave.jl") #convert weave file to markdown
makedocs()
run(`mkdocs build`)


text:
module BenchmarkConfigSweepsTests

include("test_configurations.jl")
include("test_info.jl")
include("test_examples.jl")
include("test_aqua.jl")

end  # module BenchmarkConfigSweepsTests


text:
// Declaration of multiple variables of the same type in one statement:

int main () {
  int x, y;
  x = 45;
  y = -36;
  printInt(x);
  printInt(y);
  return 0 ;

}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
module RunnerTests

using Runner
using Test

# Re-include the packages so that later scripts will not throw the Undefined
# error
@runit "common.jl"

# run simple argument
@runit "./simplesum.jl 12 13 25"
@runit "./simplesum.jl 12   23    35"  # run with extra space

# examples with ArgParse
@runit "argparsesum.jl"
@runit "argparsesum.jl --num1 33 --num2 9 --ans 42"
@runit "argparsesum.jl --num1 0 --ans 42"

end


text:
module Winds

export main

using LibSerialPort, Blink, Interact, COBS

function main(; min::Int = 0, max::Int = 255, levels = min:max, baud::Int = 9600)
    @assert min < max "minimum, $min, must be smaller than maximum, $max"
    @assert 0 ≤ min ≤ 255 "minimum must be between 0 and 255"
    @assert 0 ≤ max ≤ 255 "maximum must be between 0 and 255"

    ports = get_port_list()
    if isempty(ports)
        error("no ports were detected...")
    end
    dd = dropdown(ports)
    ok = button("OK")
    w = Window()
    body!(w, hbox("Port", dd, ok))
    on(ok) do _
        serialport = open(dd[], baud)
        int = rangepicker(levels, value = [minimum(levels)])
        on(int) do i
            encode(serialport, UInt8(i[]))
        end
        body!(w, hbox("Intensity", int))
    end
end

end # module


text:
###########################################################################
#  Adaptive Hamiltonian Monte Carlo with dual averaging (HMCDA)
#
#  Reference: Hoffman M.D, Gelman A.The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.
#  arXiv, 2011
#
#  Parameters:
#    - shrinkage: parameter that controls the amount of shrinkage towards mu
#    - t0: parameter that stabilizes the initial iterations of the algorithm
#    - step: parameter for setting step size scheduling (t^-step)
#    - len: target simulation length
#    - rate: target HMC acceptance rate
#
###########################################################################

export HMCDA

println("Loading HMCDA(rate, len, shrinkage, t0, step) sampler")

###########################################################################
#                  HMC type
###########################################################################

immutable HMCDA <: MCMCSampler
  rate::Float64
  len::Float64
  shrinkage::Float64
  t0::Float64
  step::Float64
  storeLeaps::Bool

  function HMCDA(rate::Float64, len::Float64, shrinkage::Float64, t0::Float64, step::Float64, storeLeaps::Bool)
    @assert 0. < rate < 1. "Target acceptance rate ($rate) should be between 0 and 1"
    @assert len > 0 "len parameter of HMCDA sampler ($len) must be non-negative"
    @assert shrinkage > 0. "shrinkage parameter of HMCDA sampler ($shrinkage) must be positive"
    @assert t0 >= 0 "t0 parameter of HMCDA sampler ($t0) must be non-negative"

    new(rate, len, shrinkage, t0, step, storeLeaps)
  end
end

HMCDA(;rate::Float64=0.65, len::Float64=2., shrinkage::Float64=0.05, t0::Float64=10., step::Float64=0.75, storeLeaps::Bool=false) =
  HMCDA(rate, len, shrinkage, t0, step, storeLeaps)

###########################################################################
#                  HMCDA task
###########################################################################

####  Helper functions and types for HMCDA sampling task

function initializeHMCDAStep(model::MCMCModel, sample::HMCSample)
  leapStep::Float64 = 1

  # Determine towards which direction leap step will move
  s = leapfrog(sample, leapStep, model.evalallg)
  update!(s)
  p = exp(s.H-sample.H)
  a = 2*(p>0.5)-1

  # Keep moving leap step in that direction until acceptprob crosses 0.5
  while p^a > 2^(-a)
    leapStep = leapStep*2^a
    s = leapfrog(sample, leapStep, model.evalallg)
    update!(s)
    p = exp(s.H-sample.H)
  end

  return leapStep
end

####  HMCDA task
function SamplerTask(model::MCMCModel, sampler::HMCDA, runner::MCMCRunner)
  local state0
  local dualH
  local p, i
  local nLeaps, leapStep, dualLeapStep
  local leapState0, leapStates

  @assert hasgradient(model) "HMCDA sampler requires model with gradient function"

  # hook inside Task to allow remote resetting
  task_local_storage(:reset,
    (resetPars::Vector{Float64}) -> (state0 = HMCSample(copy(resetPars)); calc!(state0, model.evalallg)))

  # initialization
  state0 = HMCSample(copy(model.init))
  calc!(state0, model.evalallg)
  @assert isfinite(state0.logTarget) "Initial values out of model support, try other values"

  state0.m = randn(model.size)
  leapStep = initializeHMCDAStep(model, state0)
  mu = log(10*leapStep)
  dualLeapStep = 1.
  dualH = 0.

  #  main loop
for i in 1:Inf
    local j, state

    state0.m = randn(model.size)
    update!(state0)
    state = deepcopy(state0)

    nLeaps = max(1, round(sampler.len/leapStep))

    if !sampler.storeLeaps
      for j = 1:nLeaps
        state = leapfrog(state, leapStep, model.evalallg)
      end
    else
      leapStates = Array(HMCSample, nLeaps+1)
      leapStates[1] = deepcopy(state0)
      for j = 2:nLeaps+1
        state = leapfrog(state, leapStep, model.evalallg)
        leapStates[j] = deepcopy(state)
      end
    end

    # accept if new is good enough
    p = min(1, exp(state0.H-state.H))
    if rand() < p
      diagnostics = (!sampler.storeLeaps ? {"accept" => true} : {"accept" => true, "leaps" => leapStates})
      ms = MCMCSample(state.pars, state.logTarget, state.grad, state0.pars, state0.logTarget, state0.grad, diagnostics)
      produce(ms)
      state0 = deepcopy(state)
    else
      diagnostics = (!sampler.storeLeaps ? {"accept" => false} : {"accept" => false, "leaps" => leapStates})
      ms = MCMCSample(state0.pars, state0.logTarget, state0.grad, state0.pars, state0.logTarget, state0.grad,
        diagnostics)
      produce(ms)
    end

    if i < runner.burnin
      eta = 1/(i+sampler.t0)
      dualH = (1-eta)*dualH+eta*(sampler.rate-p)
      leapStep = exp(mu - sqrt(i)*dualH/sampler.shrinkage)
      eta = i^(-sampler.step)
      dualLeapStep = exp((1-eta)*log(dualLeapStep)+eta*log(leapStep))
    else
      leapStep = dualLeapStep
    end
  end
end


text:
using Documenter
using CommutativeFunctions

DocMeta.setdocmeta!(CommutativeFunctions, :DocTestSetup, :(using CommutativeFunctions); recursive=true)

makedocs(;
    modules=[CommutativeFunctions],
    authors="Jishnu Bhattacharya",
    repo="https://github.com/jishnub/CommutativeFunctions.jl/blob/{commit}{path}#L{line}",
    sitename="CommutativeFunctions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jishnub.github.io/CommutativeFunctions.jl",
        assets=String[],
    ),
    pages=[
        "API" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jishnub/CommutativeFunctions.jl",
    devbranch = "main",
)


text:
@testset "independent" begin
    outdim = 3
    x = KernelFunctions.MOInputIsotopicByOutputs([rand(5) for _ in 1:4], outdim)
    y = KernelFunctions.MOInputIsotopicByOutputs([rand(5) for _ in 1:4], outdim)
    z = KernelFunctions.MOInputIsotopicByOutputs([rand(5) for _ in 1:2], outdim)

    xIF = KernelFunctions.MOInputIsotopicByFeatures(x.x, outdim)
    yIF = KernelFunctions.MOInputIsotopicByFeatures(y.x, outdim)
    zIF = KernelFunctions.MOInputIsotopicByFeatures(z.x, outdim)

    k = IndependentMOKernel(GaussianKernel())
    @test k isa IndependentMOKernel
    @test k isa MOKernel
    @test k isa Kernel
    @test k.kernel isa Kernel

    @test kernelmatrix(k, x, y) == kernelmatrix(k, collect(x), collect(y))

    ## accuracy
    KernelFunctions.TestUtils.test_interface(k, x, y, z)
    KernelFunctions.TestUtils.test_interface(k, xIF, yIF, zIF)

    # type stability (maybe move to test_interface?)
    x2 = MOInput(rand(Float32, 4), 2)
    @test k(x2[1], x2[2]) isa Float32
    @test k(x2[1], x2[1]) isa Float32
    @test eltype(typeof(kernelmatrix(k, x2))) <: Float32

    @test string(k) ==
          "Independent Multi-Output Kernel\n" *
          "\tSquared Exponential Kernel (metric = Euclidean(0.0))"
end


text:
const ∇ = ContMechTensors.gradient
const Δ = ContMechTensors.hessian

function Ψ(C, μ, Kb)
    detC = det(C)
    J = sqrt(detC)
    Ĉ = detC^(-1/3)*C
    return 1/2*(μ * (trace(Ĉ)- 3) + Kb*(J-1)^2)
end

function S(C, μ, Kb)
    I = one(C)
    J = sqrt(det(C))
    invC = inv(C)
    return μ * det(C)^(-1/3)*(I - 1/3*trace(C)*invC) + Kb*(J-1)*J*invC
end

const μ = 1e10;
const Kb = 1.66e11;
Ψ(C) = Ψ(C, μ, Kb)
S(C) = S(C, μ, Kb)

@testset "AD" begin
for dim in 1:3
    println("Testing AD for dim = $dim")

    F = one(Tensor{2,dim}) + rand(Tensor{2,dim});
    C = tdot(F);
    C2 = F' ⋅ F;

    @test 2∇(Ψ, C) ≈ S(C)
    @test 2∇(Ψ, C2) ≈ S(C2)

    b = rand(SymmetricTensor{2, dim})
    @test 2 * Δ(Ψ, C) ⊡ b ≈ ∇(S, C) ⊡ b
    @test 2 * Δ(Ψ, C2) ⊡ b ≈ ∇(S, C2) ⊡ b

    @test ∇(Ψ, C) ≈ ∇(Ψ, C2)
    @test ∇(S, C) ⊡ b ≈ ∇(S, C2) ⊡ b
    @test Δ(Ψ, C) ⊡ b ≈ Δ(Ψ, C2) ⊡ b

    for T in (Float32, Float64)
        A = rand(Tensor{2, dim, T})
        B = rand(Tensor{2, dim, T})
        A_sym = rand(SymmetricTensor{2, dim, T})
        B_sym = rand(SymmetricTensor{2, dim, T})
        v = rand(Vec{dim, T})
        II = one(Tensor{4, dim, T})
        II_sym = one(SymmetricTensor{4, dim, T})

        # Gradient of scalars
        @test ∇(norm, v) ≈ ∇(norm, v, :all)[1] ≈ v / norm(v)
        @test ∇(norm, v, :all)[2] ≈ norm(v)
        @test ∇(norm, A) ≈ ∇(norm, A, :all)[1] ≈ A / norm(A)
        @test ∇(norm, A, :all)[2] ≈ norm(A)
        @test ∇(norm, A_sym) ≈ ∇(norm, A_sym, :all)[1] ≈ A_sym / norm(A_sym)
        @test ∇(norm, A_sym, :all)[2] ≈ norm(A_sym)
        @test ∇(v -> 3*v, v) ≈ ∇(v -> 3*v, v, :all)[1] ≈ diagm(Tensor{2, dim}, 3.0)
        @test ∇(v -> 3*v, v, :all)[2] ≈ 3*v
        # https://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Derivatives_of_the_invariants_of_a_second-order_tensor
        I1, DI1 = A -> trace(A), A -> one(A)
        I2, DI2 = A -> 1/2 * (trace(A)^2 - trace(A⋅A)), A -> I1(A) * one(A) - A'
        I3, DI3 = A -> det(A), A -> det(A) * inv(A)'

        @test ∇(I1, A) ≈ ∇(I1, A, :all)[1] ≈ DI1(A)
        @test ∇(I1, A, :all)[2] ≈ I1(A)
        @test ∇(I2, A) ≈ ∇(I2, A, :all)[1] ≈ DI2(A)
        @test ∇(I2, A, :all)[2] ≈ I2(A)
        @test ∇(I3, A) ≈ ∇(I3, A, :all)[1] ≈ DI3(A)
        @test ∇(I3, A, :all)[2] ≈ I3(A)
        @test ∇(I1, A_sym) ≈ ∇(I1, A_sym, :all)[1] ≈ DI1(A_sym)
        @test ∇(I1, A_sym, :all)[2] ≈ I1(A_sym)
        @test ∇(I2, A_sym) ≈ ∇(I2, A_sym, :all)[1] ≈ DI2(A_sym)
        @test ∇(I2, A_sym, :all)[2] ≈ I2(A_sym)
        @test ∇(I3, A_sym) ≈ ∇(I3, A_sym, :all)[1] ≈ DI3(A_sym)
        @test ∇(I3, A_sym, :all)[2] ≈ I3(A_sym)

        # Gradient of second order tensors
        @test ∇(identity, A) ≈ ∇(identity, A, :all)[1] ≈ II
        @test ∇(identity, A, :all)[2] ≈ A
        @test ∇(identity, A_sym) ≈ ∇(identity, A_sym, :all)[1] ≈ II_sym
        @test ∇(identity, A_sym, :all)[2] ≈ A_sym
        @test ∇(transpose, A) ⊡ B ≈ ∇(transpose, A, :all)[1] ⊡ B ≈ B'
        @test ∇(transpose, A, :all)[2] ≈ A'
        @test ∇(transpose, A_sym) ⊡ B_sym ≈ ∇(transpose, A_sym, :all)[1] ⊡ B_sym ≈ B_sym'
        @test ∇(transpose, A_sym, :all)[2] ≈ A_sym'
        @test ∇(inv, A) ⊡ B ≈ ∇(inv, A, :all)[1] ⊡ B ≈ - inv(A) ⋅ B ⋅ inv(A)
        @test ∇(inv, A, :all)[2] ≈ inv(A)
        @test ∇(inv, A_sym) ⊡ B_sym ≈ ∇(inv, A_sym, :all)[1] ⊡ B_sym ≈ - inv(A_sym) ⋅ B_sym ⋅ inv(A_sym)
        @test ∇(inv, A_sym, :all)[2] ≈ inv(A_sym)

        # Hessians of scalars
        @test Δ(norm, A).data ≈ Δ(norm, A, :all)[1].data ≈ vec(ForwardDiff.hessian(x -> sqrt(sumabs2(x)), A.data))
        @test Δ(norm, A, :all)[2].data ≈ vec(ForwardDiff.gradient(x -> sqrt(sumabs2(x)), A.data))
        @test Δ(norm, A, :all)[3] ≈ norm(A)
    end
end
end # testset


text:
using MyExample
using Test

my_f(2,1)

@testset "MyExample.jl" begin
    my_f(2,1)
end


text:
using FastCombinations
import FastCombinations: next_states

using Base.Test

@test next_states(1:6, 1, Val{0}) == ()
@test next_states(1:6, 1, Val{1}) == (1,)
@test next_states(1:6, 1, Val{6}) == (1, 2, 3, 4, 5, 6)

@test_throws ErrorException Combinations{0}(1:1)

@test first(Combinations{1}(1:6)) == (1,)
@test first(Combinations{2}(1:6)) == (1,2)
@test first(Combinations{3}(1:6)) == (1,2,3)
@test first(Combinations{4}(1:6)) == (1,2,3,4)
@test first(Combinations{5}(1:6)) == (1,2,3,4,5)
@test first(Combinations{6}(1:6)) == (1,2,3,4,5,6)
@test collect(Combinations{7}(1:6)) == NTuple{7, Int}[]

@test first(MultiCombinations{1}(1:6)) == (1,)
@test first(MultiCombinations{2}(1:6)) == (1,1)
@test first(MultiCombinations{3}(1:6)) == (1,1,1)
@test first(MultiCombinations{4}(1:6)) == (1,1,1,1)
@test first(MultiCombinations{5}(1:6)) == (1,1,1,1,1)
@test first(MultiCombinations{6}(1:6)) == (1,1,1,1,1,1)

# TODO test more complex iterators

@test length(Combinations{1}(1:1)) == 1 == length(collect(Combinations{1}(1:1)))
@test length(Combinations{1}(1:2)) == 2 == length(collect(Combinations{1}(1:2)))
@test length(Combinations{1}(1:3)) == 3 == length(collect(Combinations{1}(1:3)))

@test length(Combinations{2}(1:2)) == 1 == length(collect(Combinations{2}(1:2)))
@test length(Combinations{2}(1:3)) == 3 == length(collect(Combinations{2}(1:3)))
@test length(Combinations{2}(1:4)) == 6 == length(collect(Combinations{2}(1:4)))

@test length(Combinations{3}(1:3)) == 1 == length(collect(Combinations{3}(1:3)))
@test length(Combinations{3}(1:4)) == 4 == length(collect(Combinations{3}(1:4)))
@test length(Combinations{3}(1:5)) == 10 == length(collect(Combinations{3}(1:5)))

@test length(MultiCombinations{1}(1:1)) == 1  == length(collect(MultiCombinations{1}(1:1)))
@test length(MultiCombinations{1}(1:2)) == 2  == length(collect(MultiCombinations{1}(1:2)))
@test length(MultiCombinations{1}(1:3)) == 3  == length(collect(MultiCombinations{1}(1:3)))

@test length(MultiCombinations{2}(1:1)) == 1  == length(collect(MultiCombinations{2}(1:1)))
@test length(MultiCombinations{2}(1:2)) == 3  == length(collect(MultiCombinations{2}(1:2)))
@test length(MultiCombinations{2}(1:3)) == 6  == length(collect(MultiCombinations{2}(1:3)))

@test length(MultiCombinations{3}(1:1)) == 1  == length(collect(MultiCombinations{3}(1:1)))
@test length(MultiCombinations{3}(1:2)) == 4  == length(collect(MultiCombinations{3}(1:2)))
@test length(MultiCombinations{3}(1:3)) == 10 == length(collect(MultiCombinations{3}(1:3)))


text:
using TriMatrices
using TriMatrices.Indexing
using TriMatrices.Testing

@testset "tri and triinv" begin
	for n in [0:100; 2 .^ 10:10:30]
		t = trinum(n)

		@test t == sum(1:n)

		@test triinv(t) == n
		@test triinv(t) == n
		@test triinv_strict(t) == n
	end

	for t in [0:100; 2 .^ 10:10:30]
		n, r = triinv_rem(t)
		@test n >= 0 && (0 <= r <= n)
		@test trinum(n) + r == t
		@test triinv(t) == n
		r == 0 || @test_throws DomainError triinv_strict(t)
	end
end

@testset "Index checking and conversion" begin
	for L in LAYOUT_TYPES_P
		layout = L()
		for n in 0:10
			data, tmat = make_test_matrix(layout, n)

			# By cartesian index
			for r in 1:n, c in 1:n
				idx = CartesianIndex(r, c)
				v = tmat[r, c]

				if check_tri_index(Bool, layout, r, c)  # In data region
					@test check_tri_index(Bool, layout, idx)
					check_tri_index(layout, r, c)
					check_tri_index(layout, idx)

					i = car2lin(layout, r, c)
					@test car2lin(layout, idx) == i
					@test car2lin_unchecked(layout, r, c) == i
					@test car2lin_unchecked(layout, idx) == i

					@test tmat[r, c] == data[i]

				else  # Not stored
					@test !check_tri_index(Bool, layout, r, c)
					@test !check_tri_index(Bool, layout, idx)
					@test_throws DomainError check_tri_index(layout, r, c)
					@test_throws DomainError check_tri_index(layout, idx)
					@test_throws DomainError car2lin(layout, r, c)
					@test_throws DomainError car2lin(layout, idx)
				end
			end

			# By data index
			for i in eachindex(data)
				r, c = lin2car(layout, i)
				@test data[i] == tmat[r, c]
				@test car2lin(layout, r, c) == i
			end
		end
	end
end

@testset "tri_indices" begin
	for L in LAYOUT_TYPES_P
		layout = L()
		for n in 0:10
			data, tmat = make_test_matrix(layout, n)
			ti = tri_indices(layout, n)

			@test length(ti) == length(data)
			@test Base.IteratorEltype(ti) == Base.IteratorEltype(typeof(ti)) == Base.HasEltype()
			@test Base.IteratorSize(ti) == Base.IteratorSize(typeof(ti)) == Base.HasLength()
			!isempty(data) && @inferred first(ti)

			i = 0
			for (r, c) in tri_indices(layout, n)
				i += 1
				@test data[i] == tmat[r, c]
			end
			@test i == length(data)
		end
	end
end


text:
module AlgebraicDynamics

using Catlab
using Catlab.Theories
using Catlab.WiringDiagrams
using Catlab.Programs

include("uwd_dynam.jl")
include("dwd_dynam.jl")
include("cpg_dynam.jl")
include("trajectories.jl")
include("TontiDiagrams.jl")

end # module