text:
using LinearAlgebra
using DynamicPolynomials
using SwitchOnSafety
using Combinatorics
using SparseArrays
using JuMP, Ipopt, MosekTools,NLopt
using SpecialFunctions

include("../src/RandomTrajectories.jl")
include("../src/AlgebraicLift.jl")
include("../src/ScenarioOpti.jl")
include("../src/ProbabilisticCertificates.jl")
include("../src/WhiteBoxAnalysis.jl")

dim = 2; numMode = 2; dimIn = 1

numScen_budget = 1000

Asave = [zeros(Float64, (dim, dim)) for i=1:numMode]
Bsave = zeros(Float64, (dim, dimIn))
ratio = 1

while true

A = [2*rand(Float64, (dim, dim)).-1 for i=1:numMode]
#A = [[1 1; 0 1], [1 0; 1 1]]
#println(A)
jsrbound = white_box_jsr(A)

#println("JSR: $jsrbound")

rhomax = 0
for i in 1:numMode
    rhoA = maximum(abs.(eigvals(A[i])))
    if rhomax < rhoA
        rhomax = rhoA
    end
end
#println("Eigenvalue max: $rhomax")

if jsrbound > 1 #rhomax < jsrbound*0.9 &&

B = 2*rand(Float64, (dim, dimIn)).-1
gaTrue,K = white_box_stabilization_quad(A,B)
#println(K)
if isempty(K)
    println("A: $A")
    println("B: $B")
end
jsrboundclosed = white_box_jsr([Ai+B*K for Ai in A])

if gaTrue/jsrboundclosed > ratio
    global Asave = A
    global Bsave = B
    global ratio = gaTrue/jsrboundclosed
    println("White-box stabilization: $gaTrue")
    println("JSR closed: $jsrboundclosed")
    println("A: $Asave")
    println("B: $Bsave")
    println("Ratio: $ratio")
end

#=
AK = [Ai+B*K for Ai in A]
rhomax_close = 0
for i in 1:numMode
    rhoA = maximum(abs.(eigvals(AK[i])))
    if rhomax_close < rhoA
        rhomax_close = rhoA
    end
end
println("Eigenvalue max closed: $rhomax_close")
=#
if jsrboundclosed < gaTrue*0.9 #jsrboundclosed > 1 &&

    println("White-box stabilization: $gaTrue")
    println("JSR closed: $jsrboundclosed")
    println("A: $A")
    println("B: $B")

#=
gamma1: 1.091098994365357 || K1: [0.4450968546375577 -0.0067944232635188 -0.07242928507366847; -1.0901373033935469 -0.6362471475237879 -0.6546310496281078]
gamma2: 1.0043051965701504 || K2: [-0.017481932047868836 -0.019322752052078965 -0.13267775176698232; -0.16816570986167362 -0.16848413512732924 -0.1473949009786114]
gamma3: 0.9319170550393475 || K3: [0.2185895965435866 0.0828529524089513 -0.10940266352022417; -0.4330159942259682 -0.2689004267861104 -0.36664805316411]
A: [[-0.035713556139526315 -0.46101936337212224 0.5758464004636243; 0.5365323737379186 -0.4238668271399051 -0.15898425119264026; 0.8419591377106763 -0.12433448826055438 0.08762970902198686], [0.7536818225043662 0.8448230972467972 -0.5723985482099347; -0.01764965739578206 -0.726360371394331 0.919109231975153; -0.7616763202127323 -0.45233767277561787 0.8890705069340319], [0.9675539515043958 0.33014694100762054 0.7845206390342243; -0.7664059406402712 -1.4599068229700047 -1.374238975340281; -0.9919933573988242 0.032139509230417085 -0.7241320070196489]]
B: [-0.9738898776106395 0.2430937821271404; -0.5594898459836832 -0.8567879128803524; 0.968034729630463 0.23321416718508736]

gamma1: 1.4771721521387928 || K1: [1.5837457868149079 -0.4143607875557178 -0.8926641693341837; 2.2458480283900912 -0.49739409437647436 0.12379794005429769]
gamma2: 1.1284125028892853 || K2: [1.169458826476936 1.0133947633814673 -0.483014379311149; 2.141282305221211 0.557871222870866 0.4928558882834438]
gamma3: 0.9586177379188866 || K3: [1.4852632106087928 0.768086923841151 -0.7357944137203406; 2.569449416459834 1.45216594818393 -0.6464578631991]
A: [[0.30711010163057084 0.35743471478281297 -0.4165467097338529; -1.2090478525084418 -1.0889753204734136 0.6606482717761146; 1.4916720938947043 0.712750594788254 -0.9494855947533352], [0.07478000105066651 1.2354131020593702 1.3717629013640442; -0.800207539428841 -1.2193833745904425 0.6737888109652435; 0.833849823841502 0.5815076492256317 0.31055833688134005], [0.43205427033385724 -1.2518407274846475 1.285786041398043; -0.8350914889424306 1.3566196051386124 0.5088045250037383; 1.0077852528903604 1.4886588615734646 1.2656260379533943]]
B: [0.8021437835254663 -0.6833724710101055; 0.6501895221291156 -0.07049181782485858; 0.2859499732121846 -0.6349752471957584]

gamma1: 1.018066677187341 || K1: [2.584474642796141 -0.24162143938208935 0.099462950616573; -1.2867193442164662 0.13807898609460512 -0.19437308433321066]
gamma2: 0.9832197222746352 || K2: [1.1463420634669867 -0.3221344185415759 -0.43056835102089497; -0.510805330235787 0.16847443653003535 0.09598776622354863]
gamma3: 0.9253726445472554 || K3: [0.4026290483392731 0.3331629002414559 -0.4362828807087034; -0.16661137825136454 -0.2450569921574039 0.12637508003535894]
A: [[-0.30800169258330756 0.9505145074280321 -0.6259941895718011; -0.09813858442579693 0.04701685427672242 0.47999535937466775; -0.5337119123653387 -1.22263351642073 1.0250969836175772], [-0.7760395797649433
0.44231108160967914 -0.7322470977044158; -1.0076425423418636 0.8040581854952427 0.36770750935111374; -1.2752392250900568 0.1577843128284635 -0.04355988747606543], [0.8510269920574687 -0.7025377021362369 0.3662889318400304; -1.2401443442753204 0.5192000830316674 -0.8306209038179504; -1.308005763363062 0.8883944328469036 -0.7044580923902992]]
B: [-0.673706431037024 -0.9942210379640635; 0.5887891646350618 0.816544078417563; 0.5846952620521035 -0.08410769825592324]

gamma1: 1.2655581866366972 || K1: [-0.7523606179678266 -0.313855132517474 -0.18146019093303392; -1.598699415329922 0.40263117664776193 0.14878090883819198]
gamma2: 0.9418593783946065 || K2: [-1.1860831186635448 0.2616062380247517 -0.8299427094969802; -1.1238095874215457 0.43745239624058097 0.016050395816995855]
gamma3: 0.8821086726589813 || K3: [-0.8968314046121041 -0.14999943896045986 -0.6103608545100961; -1.1547588693711313 -0.33652194601647034 0.49950344977932337]
A: [[1.0595043843018859 -0.6438627791708056 0.00723522344691041; 0.2583554978549165 0.06987763899547095 -0.8988079061161514; -0.289275086715185 0.2213995411678067 -1.289514635518868], [-0.07381204441025302 0.194844991033865 -0.2028082999442229; -1.450850561219815 -0.7939044396088417 -0.22192956007240272; 0.6725331583952081 -1.449399767981382 0.08236411024717127], [1.4015273027335442 -1.4414767057887359 0.16407460676712216; 1.1077681042753422 -1.4266521231371279 -0.6003075175361281; -0.1589603499259864 0.9127473107730948 -1.0781043037842029]]
B: [0.1296346994002282 0.16628258052872047; -0.8117363096670669 0.12530345308656266; -0.8660337247723389 0.6781414601703961]

gamma1: 1.100564798387133 || K1: [-0.4754988921219234 0.9686501257953164]
gamma2: 0.9659557480414829 || K2: [-0.5608615591223438 0.6780315365664531]
A: [[1.2151079530621365 -0.33531513729513973; 1.2900511328638773 -1.3708985625995862], [-0.3408550384526632 -1.4962998078843397; 0.6784775763322353 -0.5621789240840522]]
B: [1.2212436025566022; 1.2344115037249868]

gamma1: 1.237023188862944 || K1: [1.6656380348300133 -1.5018804552577114]
gamma2: 0.8912866616247893 || K2: [1.2853564948996887 -0.9904177141477473]
A: [[-0.19602867727143058 0.9066579411992075; -1.191682425218234 -0.8165764607928059], [-1.4500928127989938 0.4987795138241713; -1.4907935439044315 0.13680174640346343]]
B: [0.5051340748481175; 0.35155906853515373]

White-box stabilization: 1.2515596315140356
gamma1: 1.21804211145584 || K1: [-0.49492771379150446 0.17138754707168224]
gamma2: 1.1025805507623097 || K2: [-0.8109269218012012 -0.22648847280226433]
A: [[0.5321293666918354 -1.4093176949123007; 0.3913118072170749 1.186645538188861], [-0.6346354185840684 0.4569919176747286; 0.9401935791262526 -0.4447425946384207]]
B: [0.48123383347242177; 0.6004781313788943]

White-box stabilization: 1.1252337292756291
gamma2: 0.9670755795170983 || K2: [-1.7872748458126402e-8 -1.2461763794928268e-7]
gamma1: 1.2610314995182843 || K1: [0.4979140858258914 -0.044031696357208806]
gamma2: 0.9670755795170983 || K2: [-1.7872748458126402e-8 -1.2461763794928268e-7]
A: [[1.1 0.0; 1.1 0.0], [0.0 1.1; 0.0 -1.1]]
B: [-0.9295420348363925; -0.7283669003038016]

********************************************************************************
White-box stabilization: 1.6996993773113758
JSR closed: 1.4442237715227517
JSR closed 1: 1.3356865279982715
JSR closed 2: 1.1956673430584541
********************************************************************************
gamma1: 1.3359320418484906 || K1: [0.9602082951808819 0.28190789741712763]
gamma2: 1.1972912108912892 || K2: [0.9473315633445621 0.4586615830219181]
A: [[0.299560690718218 0.9921279903594802; -0.5828937405932029 -1.4334283024388252], [-0.5717533245191793 0.9926090534389171; -1.0977155184013296 0.881840305882942]]
B: [0.6311748433573756; 0.9416093511508103]

White-box stabilization: 1.1572917246881822
JSR closed: 1.1572684714914292
JSR closed 1: 1.2200360005616773
JSR closed 2: 1.1606814421214997
JSR closed Prob1: 1.425911090191965
JSR closed Prob2: 1.343875912403204
********************************************************************************
jsr_bound1: 1.425911090191965 || K1: [-0.3509479741590583 0.6414243973390533]
jsr_bound2: 1.343875912403204 || K2: [-0.29654684439153317 0.02940580265883222]
A: [[0.8967668894523957 0.8591470561273686; -0.9947635042056926 0.052965678226623236], [-0.9304253310105608 0.8035784417559428; 0.28078770400919195 0.12841567064867432]]
B: [-0.3150204290150107; -0.7924848456311944]

A: [[-0.9658651798953035 0.3985096296438022; 0.4475503860449943 -0.45730946692696284], [0.30027599523474535 0.41190471605843415; 0.7561317624812265 0.6943226203170085]]
B: [0.7634030900512458; 0.7950504413443844]
Ratio: 1.0587283527849853

A: [[-0.9658651798953035 0.3985096296438022; 0.4475503860449943 -0.45730946692696284], [0.30027599523474535 0.41190471605843415; 0.7561317624812265 0.6943226203170085]]
B: [0.7634030900512458; 0.7950504413443844]
Ratio: 1.0587283527849853

A: [[0.8388521952671368 -0.14857670536251177; 0.8633807216781597 0.490362061793125], [0.6013596101878895 0.6697432384693109; -0.12770751192259544 0.9150610500819121]]
B: [0.9467520554779116; 0.20399193884377942]
Ratio: 1.1025338044487405

A: [[0.5688223165988293 -0.9319527689733245; 0.1382351159850712 0.9622194103015995], [0.6128965259152119 -0.3493383970117194; -0.6086026455428555 0.7983126792861888]]
B: [0.4470277518825365; 0.3302120500143886]
Ratio: 1.1230851016694519

A: [[-0.3551288408756461 -0.7687534611139699; 0.15585791034673502 0.888273516122466], [-0.5552855753806911 0.5799457932582954; -0.46463995630929 0.8225095618349383]]
B: [-0.9283359870409225; 0.054108218616240045]
Ratio: 1.1016410213969783

A: [[-0.8568628075440414 -0.5266755234199239; -0.8272553905086331 -0.6271478953810932], [-0.8939469134475648 -0.9649588910919555; 0.2625687202655267 0.6883472259929344]]
B: [-0.5860369597028803; 0.42375953316053083]
Ratio: 1.1009106556369916

A: [[0.901174622293865 -0.08869489383944629; 0.8269569776366761 0.993314713329033], [0.7936667657206211 0.9854194372202842; -0.06356358138135798 0.7352040314488115]]
B: [0.5927169628378124; -0.24484908800651395]
Ratio: 1.1189344932733898

A: [[0.9884483907628807 -0.7068091222615016; -0.9003002195090137 -0.990951516720938], [0.8310213698716118 -0.8061813540872644; 0.5209915287826292 -0.5808532982265562]]
B: [0.660211058916834; 0.011826520536278196]
Ratio: 1.0730366283985475

A: [[-0.7196693010722481 0.7362603872486089; 0.7329814217270463 0.09983056420902336], [-0.9755522690972618 -0.8955679827431435; 0.83195928316858 0.13739458769843527]]
B: [0.17446556981338768; 0.8261853369794627]
Ratio: 1.0970509881191342

A: [[0.27255516753923503 -0.30678093945243345; 0.27558866313527597 0.8715540933346464], [0.15704285284712594 0.09023848189832373; -0.917867260246092 0.867735961626027]]
B: [0.26096169679031345; 0.2655937741300747]
Ratio: 1.0856031657385525

A: [[-0.572890453958856 -0.7448127070893222; -0.544909303104943 -0.13907253647783913], [-0.16060049191299353 0.3285618812240414; -0.8284553971503263 0.9731152201540585]]
B: [0.23900251665502337; 0.38307159354825204]
Ratio: 1.0914956526368038

A: [[0.9862857830325393 -0.9628655906762353; 0.853552383240527 -0.44706514265408703], [0.931383926233575 -0.334587106156667; -0.35022364335672274 -0.46608495831363017]]
B: [0.30163279391747944; 0.9492599267898285]
Ratio: 1.111633220536938

A:  [-0.73239707473325 -0.7759541512450001; 0.6626624924274163 0.6554974407951288]
 [0.725630573905192 0.5301052861465569; 0.5809309408385976 -0.5621245111668207]
B:  0.9411913877449489
 0.4413303767850083
Ratio: 0.9

A: [[0.20944308830065594 0.28463112922261313; -0.6326970820092153 0.918521512555647], [-0.981039355084623 -0.6123453579379987; -0.043365051966182566 -0.13560416398363495]]
B: [0.2509068941663237; -0.5668162475328358]
Ratio: 0.953655145424864

A: [[0.3992966486009579 -0.8825975300456066; 0.9466203607162913 0.8220671916996221], [-0.05865322607237777 -0.32416336552150593; -0.5565696843656855 -0.8985422238439993]]
B: [-0.5919429619683738; 0.18179898964335361]
Ratio: 0.9817647067324977

A: [[0.20944308830065594 0.28463112922261313; -0.6326970820092153 0.918521512555647], [-0.981039355084623 -0.6123453579379987; -0.043365051966182566 -0.13560416398363495]]
B: [0.2509068941663237; -0.5668162475328358]
Ratio: 0.953655145424864

A: [[0.7218345749085846 0.013504049625046477; -0.763247200146052 0.823150250105976], [0.726724028000425 0.8385770769048322; -0.7232303846748045 0.19537267020502425]]
B: [0.11672890767112509; 0.4312206258057798]
Ratio: 0.9461275906687189

A: [[-0.8350343770388022 -0.006107458495185014; 1.6179435194682048 0.08247284108420949], [-0.5991942074864163 -1.7792515399993318; 1.1279527670550538 0.9404192977577894]]
B: [0.9917434505323639; 1.3179954593019438]
Ratio: 0.8497930947416346

A: [[0.94828677083526 -0.2101659042736399; 0.19875597778197385 -0.5455798309312145], [1.947971792726129 -1.970481262099395; -0.6161968652936665 0.27245749594569]]
B: [-1.2164941788317947; 0.5824551034488925]
Ratio: 0.7595095957585386

A: [[-3.204495329413295 4.309541104423182; -2.2814617881793997 -4.847416810806258], [-0.22831362951284184 1.4040030357725914; 2.026008926537825 -2.489778192500285]]
B: [-1.5225673231273773; 0.7036197784439802]
Ratio: 0.937374520680605

A: [[4.086417919490026 2.1723124078764178; -4.49798906103267 4.779648794298643], [2.5317508259661103 -4.6146818762012165; -1.8399429997449723 -4.957575716576564]]
B: [-1.7281795380537424; 4.595561599251386]
Ratio: 0.8658787049310742

A: [[1.2847752917460404 -4.294823846678941; 0.3680548348155881 0.014732983067902161], [1.1712220013664698 -0.8818060555034961; 1.6530187731608859 -2.3732736764432127]]
B: [0.32980020816855316; 2.992543632417817]
Ratio: 0.8364931366068367

A: [[-1.6253847837843716 3.595555954547308; 1.0345191229698347 1.5593831936113656], [-0.097076612905167 -4.262435989897899; 2.575877434237759 1.586551277526027]]
B: [0.468806484246338; 2.756260070416152]
Ratio: 0.8876232465445488

A: [[-2.486281431248263 4.907579862031962; 1.5579585224383834 -3.287845396063369], [-0.2529089238709137 -1.0602793211568318; 2.3880838078084636 4.924238633856245]]
B: [1.1435267881128874; 1.947585246102994]
Ratio: 0.7533669990528284

A: [[4.241207142321246 3.813148958357946; 1.0459816320091715 2.1487703631884276], [3.333688705367571 -2.0158314250817178; 2.236859420992654 -3.5707586517512335]]
B: [0.07299618929951635; -3.2472806306887603]
Ratio: 0.8217387569574383

A: [[0.7246240890812068 -4.511926717076897; -0.28032498303302944 -3.665055499040828], [3.2394719689128273 -3.810821649087668; 2.4802074166290495 -3.535887169447751]]
B: [-0.8738131428454077; 1.7296575068801907]
Ratio: 0.7169505442861009

A: [[3.7950797471863176 0.4776991132799129; -3.2639167186095897 -3.617864308316223], [4.107251735332657 -0.9735634844167329; 1.8044512593033746 -2.127808228861514]]
B: [-2.7182973303240665; 2.368807412608545]
Ratio: 0.8283989329939093

A: [[3.1223586844560174 2.019655001919607; -2.7795881425176483 -0.7947275320165215], [-3.2014387931295274 0.8119510361930544; 4.305156309944529 -1.7624449098726247]]
B: [3.820940923858389; 3.2364324136157165]
Ratio: 0.8380238682087968

A: [[-4.208828438804413 -0.6455242985703507; 3.115315451439834 4.742156375441798], [-4.261050053027935 1.24090422552195; -4.849285702257447 -3.016047020102335]]
B: [-1.1599691949525637; -2.2075460392726143]
Ratio: 0.8047260384374686

A: [[0.9223152631016855 0.22822448026393616; -1.2757222844985416 4.580661037121752], [-2.5569964218018715 0.5024167836266642; 2.8555349168411155 0.7017623628006326]]
B: [-3.050405520291193; 2.607763730803234]
Ratio: 0.7081499423232154

A: [[-2.9396466368033702 -3.675273760507267; 1.887117369477858 2.6726150387995045], [-4.109170717351384 3.0318283002604502; 0.737155993689651 -3.7320488530313645]]
B: [0.2242947168125191; -3.6133508894644106]
Ratio: 0.7727224040648161

A: [[-1.6855693079222966 -0.16653425759574958; 0.7785077624686458 -1.632104413194857], [-0.2915062738633578 -3.2823565046218084; 3.9760565120317626 -0.022743066222108865]]
B: [0.19751363983838743; 0.8639638674096854]
Ratio: 0.8871790438060302

A: [[-0.8666413136853119 0.9409336940393072; -4.360436997796078 -0.9105280420072681], [-1.0605123751726442 -2.170579522714806; 2.7367406021194824 4.504726567417791]]
B: [-4.957818693521263; 4.910396744956589]
Ratio: 0.9279979608071133

A: [[0.4169346539518015 -0.09626692422848926; 4.265499661753648 -4.118719880693799], [-2.1543572032879177 2.4029580840965226; 0.08125960691302048 -3.470598928344031]]
B: [-4.911690113488434; 1.893905664087379]
Ratio: 0.8786184998266948

=#

(state0_budget,state_budget) = generate_trajectories(1,A,numScen_budget)

batchsize = numScen_budget

K0 = zeros(dimIn,dim)

K1,ga1,P1 = soslyap_alternating(state0_budget,state_budget;B=B,d=1,batchsize=batchsize,K0=K0,tol=1e-3)
println("gamma1: $ga1 || K1: $K1 ")
Aclose1 = [Ai+B*K1 for Ai in A]
jsrboundclose1 = white_box_jsr(Aclose1)
if abs(ga1-gaTrue) < 0.1

    #K2,jsr_bound2 = probabilistc_stability_certificate(state0_budget,state_budget;B=B,numMode=numMode,d=2,batchsize=1000,K0=K0,beta=0.1,tol=1e-3)
    K2,ga2,P2 = soslyap_alternating(state0_budget,state_budget;B=B,d=2,batchsize=batchsize,K0=K0,tol=1e-3)
    println("gamma2: $ga2 || K2: $K2")
    Aclose2 = [Ai+B*K2 for Ai in A]
    jsrboundclose2 = white_box_jsr(Aclose2)

    println(repeat('*', 80))
    println("White-box stabilization: $gaTrue")
    println("JSR closed: $jsrboundclosed")
    println("JSR closed 1: $jsrboundclose1")
    println("JSR closed 2: $jsrboundclose2")

    println("JSR closed Prob1: $ga1")
    println("JSR closed Prob2: $ga2")

    println(repeat('*', 80))

if ga2 < 0.9*ga1 #&& jsr_bound2 < 1
    println("jsr_bound1: $ga1 || K1: $K1 ")
    println("jsr_bound2: $ga2 || K2: $K2")
    println("A: $A")
    println("B: $B")
    break
end
end
end
end

end

max_stars_repo_name: zhemingwang/DataDrivenSwitchControl

max_issues_repo_name: zhemingwang/DataDrivenSwitchControl

max_forks_repo_name: zhemingwang/DataDrivenSwitchControl


text:
@testset verbose=true "Conversion" begin
    @testset verbose=true "$T" for T in FloatTypes
        using EllipsisNotation
        using Random
        Random.seed!(4321)

        @testset "Array types" begin
            dumb_from_float_array(f) = mapslices(x->Quaternion{T}(x...), f, dims=(1,))[1, ..]
            dumb_from_float_array(f::Vector) = Quaternion{T}(f...)
            function dumb_to_float_array(q)
                f = Array{T}(undef, (4, size(q)...))
                f[1, ..] = getproperty.(q, :w)
                f[2, ..] = getproperty.(q, :x)
                f[3, ..] = getproperty.(q, :y)
                f[4, ..] = getproperty.(q, :z)
                f
            end
            for dims in [(), (5,), (4, 3), (2, 3, 4)]
                if dims == ()
                    q = [randn(Quaternion{T})]
                else
                    q = randn(Quaternion{T}, dims)
                end
                f = to_float_array(q)
                @test size(f) == (4, size(q)...)
                @test f == dumb_to_float_array(q)
                @test q == dumb_from_float_array(f)
                @test q == from_float_array(f)
            end
            q = randn(Quaternion{T})
            f = to_float_array(q)
            @test f isa Vector{T}
            @test eltype(f) === T
            @test size(f) == (4,)
            @test f == q.components
        end

        @testset "Euler angles" begin
            N = 5_000

            random_angles = [2π .* rand(T, 3) .- π for _ in 1:5000]
            for i in 1:N
                α, β, γ = random_angles[i]
                q1 = from_euler_angles(α, β, γ)
                q2 = exp(imz*α/2) * exp(imy*β/2) * exp(imz*γ/2)
                @test q1 ≈ q2 atol=10eps(T)
            end
            q1 = from_euler_angles.(random_angles)
            q2 = map((αβγ)->(exp(αβγ[1]*imz/2)*exp(αβγ[2]*imy/2)*exp(αβγ[3]*imz/2)), random_angles)
            @test maximum(abs, q1 .- q2) < 10eps(T)

            random_rotors = randn(Rotor{T}, N)
            for i in 1:N
                q1 = random_rotors[i]
                q2 = from_euler_angles(to_euler_angles(q1))
                @test q1 ≈ q2 atol=250eps(T)
            end
            q1 = random_rotors
            q2 = from_euler_angles.(to_euler_angles.(random_rotors))
            @test maximum(abs, q1 .- q2) < 250eps(T)
        end

        @testset "Euler phases" begin
            N = 5_000
            ϵ = (T === Float16 ? 20eps(T) : 10eps(T))

            dumb_to_euler_phases(α, β, γ) = [exp(α*im), exp(β*im), exp(γ*im)]

            random_angles = [[3(π*rand(T))-π, π*rand(T), 3(π*rand(T))-π] for _ in 1:5000]
            for i in 1:N
                α, β, γ = random_angles[i]
                z1 = dumb_to_euler_phases(α, β, γ)

                q1 = from_euler_angles(α, β, γ)
                q2 = from_euler_phases(z1...)
                @test min(abs(q1-q2), abs(q1+q2)) < ϵ
                q3 = from_euler_phases(z1)
                @test q2 == q3

                z2 = to_euler_phases(q1)
                @test z1[2] ≈ z2[2] atol=ϵ
                if abs(β) < ϵ  # Typically only happens with Float16
                    @test z1[2] ≈ one(T) atol=ϵ
                    @test z2[2] ≈ one(T) atol=ϵ
                    @test z1[1]*z1[3] ≈ z2[1]*z2[3] atol=ϵ
                elseif abs(β-π) < ϵ  # Typically only happens with Float16
                    @test z1[2] ≈ -one(T) atol=ϵ
                    @test z2[2] ≈ -one(T) atol=ϵ
                    @test z1[1]*z1[3] ≈ z2[1]*z2[3] atol=ϵ
                else
                    @test z1[1] ≈ z2[1] atol=ϵ
                    @test z1[3] ≈ z2[3] atol=ϵ
                end

                q4 = from_euler_angles(α, zero(T), γ)
                @test iszero(q4.x)
                @test iszero(q4.y)
                z3 = dumb_to_euler_phases(α, zero(T), γ)
                z4 = to_euler_phases(q4)
                @test z3[2] ≈ one(T) atol=ϵ
                @test z4[2] ≈ one(T) atol=ϵ
                @test z3[1]*z3[3] ≈ z4[1]*z4[3] atol=ϵ

                q6 = from_euler_angles(α, T(π), γ)
                @test q6.w ≈ zero(T) atol=ϵ
                @test q6.z ≈ zero(T) atol=ϵ
                z5 = dumb_to_euler_phases(α, T(π), γ)
                z6 = to_euler_phases(q6)
                @test z5[2] ≈ -one(T) atol=ϵ
                @test z6[2] ≈ -one(T) atol=ϵ
                @test z5[1]*conj(z5[3]) ≈ z6[1]*conj(z6[3]) atol=(T === Float16 ? 20ϵ : ϵ)

                q7 = Quaternion(0, q6.x, q6.y, 0)
                z7 = to_euler_phases(q7)
                @test z7[2] == -one(T)
                @test z7[1]*conj(z7[3]) ≈ exp((α-γ)*im) atol=10eps(T)
            end
        end

        @testset "Spherical coordinates" begin
            N = 5_000

            random_angles = [2π .* rand(T, 3) .- π for _ in 1:5000]
            for i in 1:N
                α, β, γ = random_angles[i]
                q1 = from_spherical_coordinates(β, α)
                q2 = exp(imz*α/2) * exp(imy*β/2) * exp(imz*γ/2)
                @test q1*imz*inv(q1) ≈ q2*imz*inv(q2) atol=10eps(T)
            end

            random_rotors = randn(Rotor{T}, N)
            for i in 1:N
                q1 = random_rotors[i]
                q2 = from_spherical_coordinates(to_spherical_coordinates(q1))
                @test q1*imz*inv(q1) ≈ q2*imz*inv(q2) atol=420eps(T)
            end
        end

        @testset "Rotation matrices" begin
            if isbitstype(T)
                for _ in 1:5_000
                    q1 = randn(Rotor{T})
                    q2 = from_rotation_matrix(to_rotation_matrix(q1))
                    @test min(abs(q1-q2), abs(q1+q2)) < 50eps(T)
                end
            end
        end
    end
end

ext: jl


text:
# Julia wrapper for header: /usr/local/include/gdal_vrt.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

"""
    VRTCreate(int,
              int) -> VRTDatasetH
"""
function VRTCreate(arg1::Cint, arg2::Cint)
    ccall((:VRTCreate, libgdal), VRTDatasetH, (Cint, Cint), arg1, arg2)
end

"""
    VRTFlushCache(VRTDatasetH) -> void
"""
function VRTFlushCache(arg1::VRTDatasetH)
    ccall((:VRTFlushCache, libgdal), Void, (VRTDatasetH,), arg1)
end

"""
    VRTSerializeToXML(VRTDatasetH,
                      const char *) -> CPLXMLNode *
"""
function VRTSerializeToXML(arg1::VRTDatasetH, arg2)
    ccall((:VRTSerializeToXML, libgdal), Ptr{CPLXMLNode}, (VRTDatasetH, Cstring), arg1, arg2)
end

"""
    VRTAddBand(VRTDatasetH,
               GDALDataType,
               char **) -> int
"""
function VRTAddBand(arg1::VRTDatasetH, arg2::GDALDataType, arg3)
    ccall((:VRTAddBand, libgdal), Cint, (VRTDatasetH, GDALDataType, Ptr{Cstring}), arg1, arg2, arg3)
end

"""
    VRTAddSource(VRTSourcedRasterBandH,
                 VRTSourceH) -> CPLErr
"""
function VRTAddSource(arg1::VRTSourcedRasterBandH, arg2::VRTSourceH)
    ccall((:VRTAddSource, libgdal), CPLErr, (VRTSourcedRasterBandH, VRTSourceH), arg1, arg2)
end

"""
    VRTAddSimpleSource(VRTSourcedRasterBandH,
                       GDALRasterBandH,
                       int,
                       int,
                       int,
                       int,
                       int,
                       int,
                       int,
                       int,
                       const char *,
                       double) -> CPLErr
"""
function VRTAddSimpleSource(arg1::VRTSourcedRasterBandH, arg2::GDALRasterBandH, arg3::Cint, arg4::Cint, arg5::Cint, arg6::Cint, arg7::Cint, arg8::Cint, arg9::Cint, arg10::Cint, arg11, arg12::Cdouble)
    ccall((:VRTAddSimpleSource, libgdal), CPLErr, (VRTSourcedRasterBandH, GDALRasterBandH, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cstring, Cdouble), arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12)
end

"""
    VRTAddComplexSource(VRTSourcedRasterBandH,
                        GDALRasterBandH,
                        int,
                        int,
                        int,
                        int,
                        int,
                        int,
                        int,
                        int,
                        double,
                        double,
                        double) -> CPLErr
"""
function VRTAddComplexSource(arg1::VRTSourcedRasterBandH, arg2::GDALRasterBandH, arg3::Cint, arg4::Cint, arg5::Cint, arg6::Cint, arg7::Cint, arg8::Cint, arg9::Cint, arg10::Cint, arg11::Cdouble, arg12::Cdouble, arg13::Cdouble)
    ccall((:VRTAddComplexSource, libgdal), CPLErr, (VRTSourcedRasterBandH, GDALRasterBandH, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cdouble, Cdouble, Cdouble), arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12, arg13)
end

"""
    VRTAddFuncSource(VRTSourcedRasterBandH,
                     VRTImageReadFunc,
                     void *,
                     double) -> CPLErr
"""
function VRTAddFuncSource(arg1::VRTSourcedRasterBandH, arg2::VRTImageReadFunc, arg3, arg4::Cdouble)
    ccall((:VRTAddFuncSource, libgdal), CPLErr, (VRTSourcedRasterBandH, VRTImageReadFunc, Ptr{Void}, Cdouble), arg1, arg2, arg3, arg4)
end


text:
module FortranStrings

export AbstractFortranString, FortranString, REPEAT, LEN, LEN_TRIM, TRIM, LNBLNK, LSAME, INDEX, SCAN, @F_str, @F8_str

import Base.Broadcast: AbstractArrayStyle, DefaultArrayStyle, Broadcasted

abstract type AbstractFortranString{CharType} <: AbstractVector{CharType} end

"""
    FortranString{CharType}

Datatype for emulating FORTRAN `CHARACTER*N` strings.

This type is a wrapper for a `Vector` with elements of type `CharType`.
Some functions and broadcasting added for limited strings emulation.

Assignment should be done via broadcasting:
```jldoctest
julia> s = F8"abc"; s .= "AB" * F"CD"; s
F8"ABC"

julia> s = F"abc"; s[2:end] .= "ABC"; s
F"aAB"
```
Element-wise assignment is also supported, but it is not FORTRAN alike.
Element comparison
```jldoctest
julia> s = F8"abc"; s[1:1] == 'a'
true

julia> s = F8"abc"; (s[1:1] == 'a', s[2:2] == "b")
(true, true)
```

See also `@F_str` for quoting.
"""
struct FortranString{CharType} <: AbstractFortranString{CharType}
    data :: Vector{CharType}
end

const ForStr = FortranString

(::Type{ForStr{T}})(ch::AbstractChar) where {T} = ForStr{T}([T(ch)])
(::Type{ForStr{T}})(::UndefInitializer, n) where {T} = ForStr{T}(Array{T}(undef, n))
function (::Type{ForStr{T}})(s::AbstractString, flags::AbstractString = "") where {T}
    if length(flags) == 0
        return ForStr{T}(map(T, collect(s)))
    elseif lowercase(flags) == "q"
        return ForStr{T}(map(T, collect(replace(s, "''" => "'"))))
    elseif lowercase(flags) == "qq" || lowercase(flags) == "d"
        return ForStr{T}(map(T, collect(replace(s, "\"\"" => "\""))))
    else
        throw(ArgumentError("unknown FortranString flag: $flags"))
    end
end

Base.convert(A::Type{AbstractFortranString}, s::AbstractString) = A(s)
Base.convert(::Type{String}, fs::AbstractFortranString) = join(map(Char, fs))

Base.promote_rule(::Type{Union{AbstractChar, AbstractString}}, B::Type{AbstractFortranString}) = B
Base.promote_rule(A::Type{AbstractFortranString}, ::Type{Union{AbstractChar, AbstractString}}) = A

Base.string(fs::ForStr) = convert(String, fs)

Base.:*(fs1::ForStr{T1}, fs2::ForStr{T2}) where {T1} where {T2}=
    ForStr{T1}(vcat(fs1.data, map(T1, fs2.data)))
Base.:*(s::Union{AbstractChar, AbstractString}, fs::ForStr{T}) where {T} =
    ForStr{T}(vcat(ForStr{T}(s).data, fs.data))
Base.:*(fs::ForStr{T}, s::Union{AbstractChar, AbstractString}) where {T} =
    ForStr{T}(vcat(fs.data, ForStr{T}(s).data))

for o in (:(==), :<, :(<=), :>, :(>=), :(!=))
    @eval begin
        Base.$o(A::Union{AbstractChar, AbstractString}, B::ForStr) = ($o)(string(A), string(B))
        Base.$o(A::ForStr, B::Union{AbstractChar, AbstractString}) = ($o)(string(A), string(B))
        Base.$o(A::AbstractVector, B::ForStr) = ($o)(string(ForStr{Char}(A)), string(B))
        Base.$o(A::ForStr, B::AbstractVector) = ($o)(string(A), string(ForStr{Char}(B)))
        Base.$o(A::ForStr, B::ForStr) = ($o)(string(A), string(B))
    end
end

@inline Base.length(fs::ForStr) = length(fs.data)
@inline Base.ncodeunits(fs::ForStr) = length(fs)
@inline Base.codeunit(fs::ForStr{T}) where {T} = T<:Char ? UInt8 : unsigned(T)
@inline Base.sizeof(fs::ForStr) = sizeof(fs.data)
@inline Base.isvalid(fs::ForStr, i::Integer) = 1 <= i <= length(fs)

Base.IteratorSize(::Type{ForStr{T}}) where {T} = Base.HasLength()
Base.IteratorEltype(::Type{ForStr{T}}) where {T} = Base.HasEltype()

@inline Base.iterate(fs::ForStr) = iterate(fs.data)
@inline Base.iterate(fs::ForStr, state::Integer) = iterate(fs.data, state)

@inline Base.setindex!(fs::ForStr{T}, c, i::Integer) where {T} = setindex!(fs.data, T(c), i)
@inline Base.getindex(fs::ForStr{T}, i::Integer) where {T} = fs.data[i]
@inline Base.getindex(fs::ForStr{T}, r::AbstractRange{K}) where {T, K} = ForStr{T}(fs.data[r])

@inline Base.ndims(::ForStr) = 1
@inline Base.ndims(::Type{ForStr}) = 1
@inline Base.size(fs::ForStr) = size(fs.data)
@inline Base.axes(fs::ForStr) = axes(fs.data)
@inline Base.eltype(fs::ForStr{T}) where {T} = T
Base.similar(fs::ForStr{T}) where {T} = ForStr{T}(undef, length(fs))

function Base.show(io::IO, ::MIME"text/plain", fs::FortranString)
    show(io, fs)
    return
end
function Base.show(io::IO, fs::FortranString{T}) where {T}
    if T === Char
        print(io, 'F')
    elseif T === UInt8
        print(io, 'F', '8')
    else
        print(io, 'F', T)
    end
    show(io, string(fs))
    return
end
Base.print(io::IO, fs::FortranString) = print(io, string(fs))
Base.textwidth(fs::FortranString) = textwidth(string(fs))

# Broadcasting

struct FortranStringStyle <: AbstractArrayStyle{1} end

const ForStrStyle = FortranStringStyle

ForStrStyle(::Val{0}) = ForStrStyle()
ForStrStyle(::Val{1}) = ForStrStyle()
ForStrStyle(::Val{N}) where N = DefaultArrayStyle{N}()
Base.Broadcast.BroadcastStyle(::Type{<:ForStr}) = ForStrStyle()

function Base.Broadcast.instantiate(bc::Broadcasted{ForStrStyle})
    if bc.axes isa Nothing
        bcaxes = Broadcast.combine_axes(bc.args...)
    else
        bcaxes = bc.axes
        # FortranString is flexible in assignment in any direction thus any sizes are allowed
        #check_broadcast_axes(axes, bc.args...)
    end
    return Broadcasted{ForStrStyle}(bc.f, bc.args, bcaxes)
end
Base.similar(bc::Broadcasted{ForStrStyle}, ::Type{ElType}) where {ElType} = similar(bc, ElType, axes(bc))
Base.similar(::Broadcasted{ForStrStyle}, ::Type{ElType}, dims) where {ElType} = similar(ForStr{ElType}, dims)
Base.similar(::Broadcasted{ForStrStyle}, ::Type{Bool}, dims) = similar(BitArray, dims)

@inline unref(x) = isa(x, Base.RefValue{String}) ? x[] : x

const AllCharTypes = Union{AbstractChar, AbstractString, Ref{<:AbstractString}, AbstractVector{<:AbstractChar}}
const CmpOps = Union{typeof(<), typeof(<=), typeof(==), typeof(>=), typeof(>), typeof(!=)}

Base.copy(bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{Broadcasted,Broadcasted}}) =
     fortranstringbroadcast!(bc.f, similar(bc, Broadcast.combine_eltypes(bc.f, bc.args)), copy(bc.args[1]), copy(bc.args[2]))

Base.copy(bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{Broadcasted,<:AllCharTypes}}) =
     fortranstringbroadcast!(bc.f, similar(bc, Broadcast.combine_eltypes(bc.f, (bc.args[1], Char[]))), copy(bc.args[1]), unref(bc.args[2]))
Base.copy(bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{Broadcasted,<:Any}}) =
     fortranstringbroadcast!(bc.f, similar(bc, Broadcast.combine_eltypes(bc.f, bc.args)), copy(bc.args[1]), bc.args[2])

Base.copy(bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{<:AllCharTypes,Broadcasted}}) =
     fortranstringbroadcast!(bc.f, similar(bc, Broadcast.combine_eltypes(bc.f, (Char[], bc.args[2]))), unref(bc.args[1]), copy(bc.args[2]))
Base.copy(bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{<:Any,Broadcasted}}) =
     fortranstringbroadcast!(bc.f, similar(bc, Broadcast.combine_eltypes(bc.f, bc.args)), bc.args[1], copy(bc.args[2]))

Base.copy(bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{<:AllCharTypes,<:AllCharTypes}}) =
     fortranstringbroadcast!(bc.f, similar(bc, Broadcast.combine_eltypes(bc.f, (Char[], Char[]))), unref(bc.args[1]), unref(bc.args[2]))
Base.copy(bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{<:Any,<:AllCharTypes}}) =
     fortranstringbroadcast!(bc.f, similar(bc, Broadcast.combine_eltypes(bc.f, (bc.args[1], Char[]))), bc.args[1], unref(bc.args[2]))
Base.copy(bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{<:AllCharTypes,<:Any}}) =
     fortranstringbroadcast!(bc.f, similar(bc, Broadcast.combine_eltypes(bc.f, (Char[], bc.args[2]))), unref(bc.args[1]), bc.args[2])
Base.copy(bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{<:Any,<:Any}}) =
     fortranstringbroadcast!(bc.f, similar(bc, Broadcast.combine_eltypes(bc.f, bc.args)), bc.args[1], bc.args[2])

Base.copy(bc::Broadcasted{ForStrStyle, <:Any, <:CmpOps, <:Tuple{<:AllCharTypes,<:AllCharTypes}}) =
     fortranstringbroadcast!(bc.f, similar(bc, Broadcast.combine_eltypes(bc.f, (Char[], Char[]))), unref(bc.args[1]), unref(bc.args[2]))
Base.copy(bc::Broadcasted{ForStrStyle, <:Any, <:CmpOps, <:Tuple{<:Any,<:AllCharTypes}}) =
     fortranstringbroadcast!(bc.f, similar(bc, Broadcast.combine_eltypes(bc.f, (Char[], Char[]))), bc.args[1], unref(bc.args[2]))
Base.copy(bc::Broadcasted{ForStrStyle, <:Any, <:CmpOps, <:Tuple{<:AllCharTypes,<:Any}}) =
     fortranstringbroadcast!(bc.f, similar(bc, Broadcast.combine_eltypes(bc.f, (Char[], Char[]))), unref(bc.args[1]), bc.args[2])
Base.copy(bc::Broadcasted{ForStrStyle, <:Any, <:CmpOps, <:Tuple{<:Any,<:Any}}) =
     fortranstringbroadcast!(bc.f, similar(bc, Broadcast.combine_eltypes(bc.f, bc.args)), bc.args[1], bc.args[2])

Base.copyto!(dest::AbstractVector, bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{Broadcasted,Broadcasted}}) =
    # TODO: an unwanted copy is making here
    fortranstringbroadcast!(bc.f, dest, copyto!(dest, bc.args[1]), copy(bc.args[2]))

Base.copyto!(dest::AbstractVector, bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{Broadcasted,<:Any}}) =
    fortranstringbroadcast!(bc.f, dest, copyto!(dest, bc.args[1]), bc.args[2])

Base.copyto!(dest::AbstractVector, bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{<:Any,Broadcasted}}) =
    fortranstringbroadcast!(bc.f, dest, bc.args[1], copyto!(dest, bc.args[2]))

Base.copyto!(dest::AbstractVector, bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{<:Any,<:Any}}) =
    fortranstringbroadcast!(bc.f, dest, bc.args[1], bc.args[2])

Base.copyto!(dest::AbstractVector, bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{Broadcasted}}) =
    fortranstringcopyto!(bc.f, dest, copyto!(dest, first(bc.args)))

Base.copyto!(dest::AbstractVector, bc::Broadcasted{ForStrStyle, <:Any, <:Any, <:Tuple{<:Any}}) =
    fortranstringcopyto!(bc.f, dest, bc.args[1])

# special case for `String` because it has `DefaultArrayStyle{0}`
# see also broadcast.jl:852 -- `copy(bc::Broadcasted{<:AbstractArrayStyle{0}}) = bc[CartesianIndex()]`
Base.copyto!(dest::Union{SubArray{T,N,P},ForStr{T}}, bc::Broadcasted{DefaultArrayStyle{0}}) where {T,N,P<:ForStr} =
    fortranstringcopyto!(bc.f, dest, bc.args[1])

for (Tc,oc) in [(Any, :UInt64), (CmpOps, :Char)]
    for (T1,o1) in [(Any, :identity), (Ref{<:AbstractString}, :getindex)]
        for (T2,o2) in [(Any, :identity), (Ref{<:AbstractString}, :getindex)]
            """
            Apply `op` element-wise on `src1` and `src2` converting to `Char`
            in of comparison operations, and store the result in `dest`.
            """
            @eval function fortranstringbroadcast!(op::$Tc, dest, src1::$T1, src2::$T2)
                i = firstindex(dest) - 1
                l1 = length($o1(src1))
                l2 = length($o2(src2))
                if l1 == l2 == 1
                    v1 = first($o1(src1))
                    v2 = first($o2(src2))
                    dest[i+=1] = op($oc(v1), $oc(v2))
                elseif l1 == 1
                    v1 = first($o1(src1))
                    for v2 in $o2(src2)
                        dest[i+=1] = op($oc(v1), $oc(v2))
                    end
                elseif l2 == 1
                    v2 = first($o2(src2))
                    for v1 in $o1(src1)
                        dest[i+=1] = op($oc(v1), $oc(v2))
                    end
                else
                    for (v1,v2) = zip($o1(src1), $o2(src2))
                        dest[i+=1] = op($oc(v1), $oc(v2))
                    end
                end
                return dest
            end
        end
    end
end

for (T1,o1) in [(Any, :identity), (Ref{<:AbstractString}, :getindex)]
    """
    Apply `op` element-wise on `src` converted to `Char` and set the result into `dest`.
    The result vector will be stretched or shrinked to satisfy size of `dest`.
    """
    @eval function fortranstringcopyto!(op, dest, src::$T1)
        (lendest = length(dest)) == 0 && return dest
        T = eltype(dest)
        i = firstindex(dest)
        lastind = lastindex(dest)
        for v in $o1(src)
            i <= lastind || break
            dest[i] = T(op(Char(v)))
            i += 1
        end
        while i <= lastind
            dest[i] = T(' ')
            i += 1
        end
        return dest
    end
end

# Some `String` and FORTRAN functions

Base.strip(fs::AbstractFortranString) = lstrip(rstrip(fs))

Base.lstrip(fs::AbstractFortranString) =
    typeof(fs)(fs.data[something(findfirst(c->!isspace(Char(c)), fs),end+1):end])

Base.rstrip(fs::AbstractFortranString) =
    typeof(fs)(fs.data[1:something(findlast(c->!isspace(Char(c)), fs),0)])

function Base.repeat(fs::AbstractFortranString, n::Integer)
    l = length(fs)
    r = typeof(fs)(undef, n*l)
    for i = 0:n-1
        @view(r[1 + i*l:l + i*l]) .= fs
    end
    return r
end

Base.uppercase(fs::ForStr) = typeof(fs)(map(uppercase, string(fs)))
Base.lowercase(fs::ForStr) = typeof(fs)(map(lowercase, string(fs)))
Base.uppercase(fs::AbstractFortranString) = typeof(fs)(map(uppercase, Char.(fs.data)))
Base.lowercase(fs::AbstractFortranString) = typeof(fs)(map(lowercase, Char.(fs.data)))

Base.occursin(r, fs::AbstractFortranString) = occursin(r, string(fs))
Base.occursin(r::AbstractFortranString, fs::AbstractFortranString) = occursin(string(r), string(fs))
Base.occursin(r::AbstractFortranString, fs) = occursin(string(r), fs)
Base.occursin(r::Regex, fs::AbstractFortranString; offset::Integer=0) = occursin(r, string(fs), offset=offset)

Base.match(r, fs::AbstractFortranString) = match(r, string(fs))

Base.split(fs::AbstractFortranString, splitter::AbstractFortranString; limit::Integer=0, keepempty::Bool=true) =
    map(typeof(fs), split(string(fs), string(splitter), limit=limit, keepempty=keepempty))
Base.split(fs::AbstractFortranString, splitter; limit::Integer=0, keepempty::Bool=true) =
    map(typeof(fs), split(string(fs), splitter, limit=limit, keepempty=keepempty))
Base.split(fs, splitter::AbstractFortranString; limit::Integer=0, keepempty::Bool=true) =
    map(typeof(fs), split(fs, string(splitter), limit=limit, keepempty=keepempty))

"""
    REPEAT(STRING, N)

Repeat `STRING` `N` times.
"""
REPEAT(s::Union{AbstractChar,AbstractString,AbstractFortranString}, n::Integer) = repeat(s, n)

"""
    LEN(STRING)

FORTRAN compatible name alias for `length()`.
"""
LEN(s::Union{AbstractChar,AbstractString,AbstractFortranString}) = length(s)

"""
    LEN_TRIM(STRING)

Return the length of the `STRING` with trailing spaces ignored.
"""
LEN_TRIM(s::Union{AbstractChar,AbstractString,AbstractFortranString}) =
    something( findlast(c->!isspace(Char(c)), string(s)), 0 )

"""
    LNBLNK(STRING)

Alias for LEN_TRIM.
"""
LNBLNK(s::Union{AbstractString,AbstractFortranString}) = LEN_TRIM(s)

"""
    TRIM(STRING)

Return a truncated copy of the `STRING` where all trailing
spaces are removed.
"""
TRIM(s::Union{AbstractString,AbstractFortranString}) = rstrip(s)

"""
    LSAME(s1, s2)

Standard BLAS function to compare first character of two strings,
returns `bool`.
"""
LSAME(s1, s2) = uppercase(Char(first(s1))) == uppercase(Char(first(s2)))

"""
    INDEX(STRING, SUBSTRING, BACK=false)

Return the position of the start of the first occurrence of string `SUBSTRING`
as a substring in `STRING`, counting from one.

If `SUBSTRING` is not present in `STRING`, zero is returned.
If the `BACK` argument is present and true, the return value is the start of
the last occurrence rather than the first.
"""
INDEX(s::Union{AbstractString,AbstractFortranString}, substr, back::Bool) =
    first(something( back ? findlast(string(substr), string(s)) : findfirst(string(substr), string(s)), 0 ))
@inline INDEX(s::Union{AbstractString,AbstractFortranString}, substr; back::Bool=false) =
    INDEX(s, substr, back)

"""
    SCAN(STRING, CHARSET, BACK=false)

Return the position of the start of the first occurrence of string `CHARSET`
as a `CHARSET` in `STRING`, counting from one.

If `CHARSET` is not present in `STRING`, zero is returned. If the `BACK` argument
is present and true, the return value is the start of the last occurrence rather than
the first.
"""
function SCAN(s::Union{AbstractString,AbstractFortranString}, charset, back::Bool)
    if back
        i = 0
        for c in charset
            i = max( something(findlast(==(Char(c)), string(s)), i), i)
        end
        return i
    else
        i = length(s) + 1
        for c in charset
            i = min( something(findfirst(==(Char(c)), string(s)), i), i)
        end
        return i > length(s) ? 0 : i
    end
end
@inline SCAN(s::Union{AbstractString,AbstractFortranString}, charset; back::Bool=false) =
    SCAN(s, charset, back)

"""
    F"some string"

String macro for `FortranString{Char}("some string")`.

The `'d'` and `"qq"` flag can be used to emulate behaviour fortran strings in double quotes,
where the double quotes inside such strings must be doubled:

    F8"Use \\"\\"doubled\\"\\""qq

The `'q'` flag indicates the emulation of Fortran strings in single quotes:

    F8"Can''t be with only single quote inside"q

"""
macro F_str(s, flags...) ; FortranString{Char}(s, flags...) ; end

"""
    F8"some string"

String macro for compatible with FORTRAN strings `FortranString{UInt8}`.

The `'d'` and `"qq"` flag can be used to emulate behaviour fortran strings in double quotes,
where the double quotes inside such strings must be doubled:

    F8"Use \\"\\"doubled\\"\\""qq

The `'q'` flag indicates the emulation of Fortran strings in single quotes:

    F8"Can''t be with only single quote inside"q

"""
macro F8_str(s, flags...) ; FortranString{UInt8}(s, flags...) ; end

end # of module FortranStrings


text:
# Basic constructors of descriptor systems
"""
    sys = dss(A, E, B, C, D; Ts = 0, check_reg = false,
              atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ )

Create for `Ts = 0` a descriptor system model `sys = (A-λE,B,C,D)` for a continuous-time state space system of the form

    Edx(t)/dt = Ax(t) + Bu(t) ,
    y(t)      = Cx(t) + Du(t) ,

where `x(t)`, `u(t)` and `y(t)` are the system state vector, system input vector and system output vector, respectively,
for the continuous time variable `t`.

For a nonzero positive sampling time `Ts = ΔT`, the descriptor system model specifies a discrete-time state space system of the form

    Ex(t+ΔT) = Ax(t) + Bu(t)
    y(t)     = Cx(t) + Du(t)

for the discrete values of the time variable `t = 0, ΔT, 2ΔT, ...`.
Use `Ts = -1` if the sampling time is not specified. In this case, by convention  `ΔT = 1`.

For a system with zero feedthrough matrix `D`, it is possible to set `D = 0` (the scalar zero).

For a standard state space system, `E` is the identity matrix. In this case, it is possible to set `E = I` (the boolean uniform scaling).
Alternatively, use

    sys = dss(A, B, C, D; Ts = 0)

to create a standard system.

For a system corresponding to a static gain `D`, use

    sys = dss(D; Ts = 0)

It is possible to specify a descriptor system via all or part of its matrices using the form

    sys = dss(A = mat1, E = mat2, B = mat3, C = mat4, D = mat5; Ts = 0, check_reg = false,
              atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ)

where `A`, `E`, `B`, `C`, and `D` are keyword parameters set to appropriate matrix values
`mat1`, `mat2`, `mat3`, `mat4`, and `mat5`, respectively. If
some of the system matrices are omited, then zero matrices of appropriate sizes are employed instead.

It is assumed that the pencil `A-λE` is regular (i.e., `det(A-λE) ̸≡ 0`), and therefore, in the interest of efficiency,
the regularity of `A-λE` is by default _not_ tested. If `check_reg = true`, the regularity of `A-λE` is
additionally checked. In this case, the keyword arguments `atol1`, `atol2` and `rtol` specify, respectively,
the absolute tolerance for the nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A` and `E`.
The default relative tolerance is `n*ϵ`, where `n` is the order of the square matrices `A` and `E`, and  `ϵ` is the working machine epsilon.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.
"""
function dss(A::AbstractNumOrArray, E::Union{AbstractNumOrArray,UniformScaling}, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray;
             Ts::Real = 0, check_reg::Bool = false,
             atol::Real = zero(real(eltype(A))), atol1::Real = atol, atol2::Real = atol,
             rtol::Real = (typeof(A) <: Number ? 1 : min(size(A)...))*eps(real(float(one(eltype(A)))))*iszero(min(atol1,atol2)))

    T = promote_type(eltype(A), E == I ? Bool : eltype(E), eltype(B), eltype(C), eltype(D))
    check_reg && E != I && !MatrixPencils.isregular(to_matrix(T,A), to_matrix(T,E), atol1 = atol1, atol2 = atol2, rtol = rtol ) &&
                            error("The pencil A-λE is not regular")
    p = typeof(C) <: Union{Vector,Number} ? (size(A,1) <= 1 ? size(C,1) : 1) : size(C,1)
    m = typeof(B) <: Union{Vector,Number} ? 1 : size(B,2)
    return DescriptorStateSpace{T}(to_matrix(T,A), E == I ? I : to_matrix(T,E), to_matrix(T,B), to_matrix(T,C,p <=1),
                                   typeof(D) <: Number && iszero(D) ? zeros(T,p,m) : p <=1 ? to_matrix(T,D,m > p) : to_matrix(T,D), Ts)
end
function dss(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray; Ts::Real = 0)
    T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D))
    p = typeof(C) <: Union{Vector,Number} ? (size(A,1) <= 1 ? size(C,1) : 1) : size(C,1)
    m = size(B,2)
    return DescriptorStateSpace{T}(to_matrix(T,A), I, to_matrix(T,B), to_matrix(T,C,p <=1),
                                   typeof(D) <: Number && iszero(D) ? zeros(T,p,m) : p <=1 ? to_matrix(T,D,m > p) : to_matrix(T,D), Ts)
end
function dss(D::AbstractNumOrArray; Ts::Real = 0)
    D == [] && (T = Int64; return DescriptorStateSpace{T}(zeros(T,0,0),I,zeros(T,0,0),zeros(T,0,0),zeros(T,0,0),Ts))
    p, m = size(D,1),size(D,2)
    T = eltype(D)
    return DescriptorStateSpace{T}(zeros(T,0,0),I,zeros(T,0,m),zeros(T,p,0),to_matrix(T,D),Ts)
end
function dss(;A::Union{AbstractNumOrArray} = zeros(Bool,0,0), E::Union{AbstractNumOrArray,UniformScaling} = I, B::Union{AbstractNumOrArray,Missing} = missing,
              C::Union{AbstractNumOrArray,Missing} = missing, D::Union{AbstractNumOrArray,Missing} = missing,
              Ts::Real = 0, check_reg::Bool = false,
              atol::Real = zero(real(eltype(A))), atol1::Real = atol, atol2::Real = atol,
              rtol::Real = (typeof(A) <: Number ? 1 : min(size(A)...))*eps(real(float(one(eltype(A)))))*iszero(min(atol1,atol2)))
    T = Bool
    T = promote_type(T, ismissing(A) ? T : eltype(A), E == I ? Bool : eltype(E), ismissing(B) ? T : eltype(B),
                     ismissing(C) ? T : eltype(C), ismissing(D) ? T : eltype(D))
    A = to_matrix(T,A)
    n = size(A,1)
    E == I || (E = to_matrix(T,E))
    check_reg && E != I && !isregular(A, E, atol1 = atol1, atol2 = atol2, rtol = rtol ) &&
                            error("The pencil A-λE is not regular")

    ismissing(D) || (D1 = to_matrix(T,D))
    if ismissing(B)
       ismissing(D) ? m = 0 : m = size(D,2)
       B = zeros(T,n,m)
    else
       B = to_matrix(T,B)
       m = size(B,2)
    end
    if ismissing(C)
        ismissing(D) ? p = 0 : p = size(D,1)
        C = zeros(T,p,n)
    else
        C = to_matrix(T,C)
        p = size(C,1)
    end
    ismissing(D) && (D1 = zeros(T,p,m))
    return DescriptorStateSpace{T}(A, E == I ? I : E, B, C, D1, Ts)
end
"""
     sys = dss(A, E, B, F, C, G, D, H; compacted = false,
               atol = 0, atol1 = atol, atol2 = atol, atol3 = atol, rtol = min(atol1,atol2,atol3)>0 ? 0 : n*ϵ)

Construct an input-output equivalent descriptor system representation `sys = (Ad-λdE,Bd,Cd,Dd)` to a pencil based linearization
`(A-λE,B-λF,C-λG,D-λH)` satisfying

                -1                        -1
     Cd*(λEd-Ad)  *Bd + Dd = (C-λG)*(λE-A)  *(B-λF) + D-λH .

If `compacted = true`, a compacted descriptor system realization is determined by exploiting possible rank defficiencies of the
matrices `F`, `G`, and `H`. Any of the matrices `F`, `G`, and `H` can be set to `missing`.

The keyword arguments `atol1`, `atol2`, `atol3`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `F`, the absolute tolerance for the nonzero elements of `G`,
the absolute tolerance for the nonzero elements of `H`  and the relative tolerance
for the nonzero elements of `F`, `G` and `H`. The default relative tolerance is `n*ϵ`, where `n` is the size of
of `A`, and `ϵ` is the machine epsilon of the element type of `A`.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol`, `atol2 = atol` and `atol3 = atol`.
"""
function dss(A::Matrix, E::Union{Matrix,UniformScaling}, B::Matrix, F::Union{Matrix,Missing},
             C::Matrix, G::Union{Matrix,Missing}, D::Matrix, H::Union{Matrix,Missing};
             Ts::Real = 0, compacted::Bool = false,
             atol::Real = zero(real(eltype(A))), atol1::Real = atol, atol2::Real = atol, atol3::Real = atol,
             rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2,atol3)))

    T = promote_type(eltype(A), E == I ? Bool : eltype(E), eltype(B), eltype(C), eltype(D))
    T = promote_type(T, ismissing(F) ? T : eltype(F), ismissing(G) ? T : eltype(G), ismissing(H) ? T :  eltype(H))

    return DescriptorStateSpace{T}(lps2ls(A, E, B, F, C, G, D, H; compacted = compacted, atol1 = atol1, atol2 = atol2, atol3 = atol2, rtol = rtol)..., Ts)

end
# function dss(sys::PencilStateSpace; compacted::Bool = false,
#              atol1::Real = zero(real(eltype(sys.A))), atol2::Real = zero(real(eltype(sys.A))),
#              rtol::Real = (min(size(sys.A)...)*eps(real(float(one(eltype(sys.A))))))*iszero(min(atol1,atol2)))
#     return DescriptorStateSpace{eltype(sys.A)}(lps2ls(sys.A,sys.E,sys.B,sys.F,sys.C,sys.G,sys.D,sys.H; compacted = compacted, atol1 = atol1, atol2 = atol2, rtol = rtol)..., sys.Ts)
# end
# function dss(D::VecOrMat, H::VecOrMat; Ts::Real = 0)
#     typeof(D) <: Vector ? (p, m) = (length(D), 1) : (p,m) = size(D)
#     T = promote_type(eltype(D), eltype(H))
#     D1 = reshape(copy_oftype(D,T),p,m)
#     iszero(H) &&  (return DescriptorStateSpace{T}(zeros(T,0,0),I,zeros(T,0,m),zeros(T,p,0),D1,Ts))
#     if m <= p
#         A = Matrix{T}(I,2*m,2*m)
#         E = [zeros(T,m,m) I; zeros(T,m,2*m)]
#         B = [zeros(T,m,m); I]
#         C = [H zeros(T,p,m)]
#      else
#         A = Matrix{T}(I,2*p,2*p)
#         E = [zeros(T,p,p) I; zeros(T,p,2*p)]
#         B = [zeros(T,p,m); -H]
#         C = [-I zeros(T,p,p)]
#      end
#      return DescriptorStateSpace{T}(A,E,B,C,D1,Ts)
# end
"""
    sys = dss(NUM, DEN; contr = false, obs = false, noseig = false, minimal = false, fast = true, atol = 0, rtol)

Convert the rational matrix `R(λ) = NUM(λ) ./ DEN(λ)` to a descriptor system representation `sys = (A-λE,B,C,D)` such that
the transfer function matrix of `sys` is `R(λ)`.

`NUM(λ)` is a polynomial matrix of the form `NUM(λ) = N_1 + λ N_2 + ... + λ**k N_(k+1)`, for which
the coefficient matrices `N_i`, `i = 1, ..., k+1` are stored in the 3-dimensional matrix `NUM`,
where `NUM[:,:,i]` contains the `i`-th coefficient matrix `N_i` (multiplying `λ**(i-1)`).

`DEN(λ)` is a polynomial matrix of the form `DEN(λ) = D_1 + λ D_2 + ... + λ**l D_(l+1)`, for which
the coefficient matrices `D_i`, `i = 1, ..., l+1`, are stored in the 3-dimensional matrix `DEN`,
where `DEN[:,:,i]` contain the `i`-th coefficient matrix `D_i` (multiplying `λ**(i-1)`).

Alternatively, `NUM(λ)` and `DEN(λ)` can be specified as matrices of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

If `n` is the order of `A-λE`, then the computed linearization satisfies:

(1) `A-λE` is regular and `R(λ) = C*inv(λE-A)*B+D`;

(2) `rank[B A-λE] = n` (controllability) if `minimal = true` or `contr = true`;

(3) `rank[A-λE; C] = n` (observability)  if `minimal = true` or `obs = true`;

(4) `A-λE` has no non-dynamic modes if `minimal = true` or `noseig = true`.

If conditions (1)-(4) are satisfied, the realization is called `minimal` and the resulting order `n`
is the least achievable order. If conditions (1)-(3) are satisfied, the realization is called `irreducible`
and the resulting order `n` is the least achievable order using orthogonal similarity transformations.
An irreducible realization preserves the pole-zero and singular structures of `R(λ)`.

The descriptor system based realization is built using the methods described in [1] in conjunction with
pencil manipulation algorithms [2] and [3] to compute reduced order realization. These algorithms
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition, if `fast = false`.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The keyword arguments `atol` and `rtol`, specify the absolute and relative tolerances, respectively, for the
nonzero coefficients of `NUM(λ)` and `DEN(λ)`.

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).

[2] P. Van Dooreen, The generalized eigenstructure problem in linear system theory,
IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981.

[3] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
function dss(NUM::Array{T1,3},DEN::Array{T2,3};
    Ts::Real = 0, minimal::Bool = false, contr::Bool = false, obs::Bool = false, noseig::Bool = false,
    atol::Real = zero(real(T1)), rtol::Real = (min(size(NUM)...)*eps(real(float(one(T1)))))*iszero(atol)) where {T1,T2}
    sysdata = rm2ls(NUM, DEN, contr = contr, obs = obs, noseig = noseig, minimal = minimal, atol = atol, rtol = rtol)
    return DescriptorStateSpace{eltype(sysdata[1])}(sysdata[1:5]...,Ts)
end
dss(NUM::Union{VecOrMat{Polynomial{T1,X}},Polynomial{T1,X},Number,VecOrMat},
    DEN::Union{VecOrMat{Polynomial{T2,X}},Polynomial{T2,X},Number,VecOrMat}; kwargs...) where {T1,T2,X} =
    dss(poly2pm(NUM),poly2pm(DEN); kwargs...)
"""
    sys = dss(R; Ts=missing, contr = false, obs = false, noseig = false, minimal = false, fast = true, atol = 0, rtol)

Convert the rational transfer function matrix `R(λ)` to a descriptor system representation `sys = (A-λE,B,C,D)` such that
the transfer function matrix of `sys` is `R(λ)`. The resulting `sys` is a continuous-time system if `Ts = 0` or
discrete-time system if `Ts = -1` or `Ts > 0`.
If `Ts = missing`, the sampling time of `sys` is inherited from the sampling time `TRs` of the elements of `R`, unless `TRs = nothing`,
in which case `Ts = 0` is used (by default).

`R(λ)` is a matrix with rational transfer function entries (see [`RationalTransferFunction`](@ref) ) corresponding to a
multiple-input multiple-outputs system or a rational transfer function corresponding to a
single-input single-output system. The numerators and denominators of the elements of R are of type
 `Polynomial` as provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

If `n` is the order of `A-λE`, then the computed realization satisfies:

(1) `A-λE` is regular and `R(λ) = C*inv(λE-A)*B+D`;

(2) `rank[B A-λE] = n` (controllability) if `minimal = true` or `contr = true`;

(3) `rank[A-λE; C] = n` (observability)  if `minimal = true` or `contr = true`;

(4) `A-λE` has no non-dynamic modes if `minimal = true` or `noseig = true`.

If conditions (1)-(4) are satisfied, the realization is called `minimal` and the resulting order `n`
is the least achievable order. If conditions (1)-(3) are satisfied, the realization is called `irreducible`
and the resulting order `n` is the least achievable order using orthogonal similarity transformations.
An irreducible realization preserves the pole-zero and singular structures of `R(λ)`.

The underlying pencil manipulation algorithms [1] and [2] to compute reduced order realizations
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition, if `fast = false`.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The keyword arguments `atol` and `rtol`, specify the absolute and relative tolerances for the
nonzero coefficients of `R(λ)`.

[1] P. Van Dooreen, The generalized eigenstructure problem in linear system theory,
IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981.

[2] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
# function dss(R::Union{VecOrMat{<:RationalTransferFunction},RationalTransferFunction}; Ts::Union{Real,Missing} = missing,
#              minimal::Bool = false, contr::Bool = false, obs::Bool = false, noseig::Bool = false,
#              atol::Real = zero(float(real(_eltype(R)))), rtol::Real = 100*eps(one(atol))*iszero(atol))
#     ismissing(Ts) && (eltype(R) <: RationalTransferFunction ? Ts = R[1].Ts : Ts = R.Ts )
#     isnothing(Ts) && (Ts = 0)
#     dss(poly2pm(numpoly.(R)), poly2pm(denpoly.(R)); Ts = Ts, contr = contr, obs = obs,
#             noseig = noseig, minimal = minimal, atol = atol, rtol = rtol)
# end
# function dss(R::VecOrMat{RationalTransferFunction{T,X}}; Ts::Union{Real,Missing} = missing,
#              minimal::Bool = false, contr::Bool = false, obs::Bool = false, noseig::Bool = false,
#              atol::Real = zero(float(real(T))), rtol::Real = 100*eps(one(atol))*iszero(atol)) where T where X
#     length(R) == 0 && (return dss(zeros(_eltype(R),size(R)...),Ts = Ts))
#     ismissing(Ts) && (Ts = R[1].Ts)
#     isnothing(Ts) && (Ts = 0)
#     dss(poly2pm(numpoly.(R)), poly2pm(denpoly.(R)); Ts = Ts, contr = contr, obs = obs,
#             noseig = noseig, minimal = minimal, atol = atol, rtol = rtol)
# end
function dss(R::VecOrMat{<:RationalTransferFunction}; Ts::Union{Real,Missing,Nothing} = missing,
             minimal::Bool = false, contr::Bool = false, obs::Bool = false, noseig::Bool = false,
             atol::Real = zero(float(real(_eltype(R)))), rtol::Real = 100*eps(one(atol))*iszero(atol))
    length(R) == 0 && (return dss(zeros(_eltype(R),size(R)...),Ts = Ts))
    ismissing(Ts) && (Ts = R[1].Ts)
    isnothing(Ts) && (Ts = 0)
    dss(poly2pm(numpoly.(R)), poly2pm(denpoly.(R)); Ts = Ts, contr = contr, obs = obs,
            noseig = noseig, minimal = minimal, atol = atol, rtol = rtol)
end
function dss(R::VecOrMat{<:RationalFunction}; Ts::Real = 0,
             minimal::Bool = false, contr::Bool = false, obs::Bool = false, noseig::Bool = false,
             atol::Real = zero(float(real(eltype(eltype(eltype(R)))))), rtol::Real = 100*eps(one(atol))*iszero(atol))
    length(R) == 0 && (return dss(zeros(eltype(eltype(eltype(G))),size(R)...),Ts = Ts))
    dss(poly2pm(numerator.(R)), poly2pm(denominator.(R)); Ts = Ts, contr = contr, obs = obs,
            noseig = noseig, minimal = minimal, atol = atol, rtol = rtol)
end
function dss(R::RationalTransferFunction; Ts::Union{Real,Missing,Nothing} = missing,
             minimal::Bool = false, contr::Bool = false, obs::Bool = false, noseig::Bool = false,
             atol::Real = zero(float(real(_eltype(R)))), rtol::Real = 100*eps(one(atol))*iszero(atol))
    ismissing(Ts) && (Ts = R.Ts)
    isnothing(Ts) && (Ts = 0)
    dss(poly2pm(numpoly(R)), poly2pm(denpoly(R)); Ts = Ts, contr = contr, obs = obs,
            noseig = noseig, minimal = minimal, atol = atol, rtol = rtol)
end
function dss(R::RationalFunction; Ts::Real = 0,
             minimal::Bool = false, contr::Bool = false, obs::Bool = false, noseig::Bool = false,
             atol::Real = zero(float(real(eltype(eltype(R))))), rtol::Real = 100*eps(one(atol))*iszero(atol))
    dss(poly2pm(numerator(R)), poly2pm(denominator(R)); Ts = Ts, contr = contr, obs = obs,
            noseig = noseig, minimal = minimal, atol = atol, rtol = rtol)
end
"""
    sys = dss(P; Ts = 0, contr = false, obs = false, noseig = false, minimal = false, fast = true, atol = 0, rtol)

Convert the polynomial matrix `P(λ)` to a descriptor system representation `sys = (A-λE,B,C,D)` such that
the transfer function matrix of `sys` is `P(λ)`. The resulting `sys` is a continuous-time system if `Ts = 0` or
discrete-time system if `Ts = -1` or `Ts > 0`.

`P(λ)` can be specified as a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

If `d` is the degree of `P(λ)` and `n` is the order of `A-λE`, then the computed realization satisfies:

(1) `A-λE` is regular and `P(λ) = C*inv(λE-A)*B+D`;

(2) `rank[B A-λE] = n` (controllability) if `minimal = true` or `contr = true`;

(3) `rank[A-λE; C] = n` (observability)  if `minimal = true` or `contr = true`;

(4) `A-λE` has no non-dynamic modes if `minimal = true` or `noseig = true`.

If conditions (1)-(4) are satisfied, the realization is called `minimal` and the resulting order `n`
is the least achievable order. If conditions (1)-(3) are satisfied, the realization is called `irreducible`
and the resulting order `n` is the least achievable order using orthogonal similarity transformations.
An irreducible realization preserves the pole-zero and singular structures of `P(λ)`.

The underlying pencil manipulation algorithms [1] and [2] to compute reduced order realizations
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition, if `fast = false`.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The keyword arguments `atol` and `rtol`, specify the absolute and relative tolerances for the
nonzero coefficients of `P(λ)`, respectively.

[1] P. Van Dooreen, The generalized eigenstructure problem in linear system theory,
IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981.

[2] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
function dss(P::Array{T,3};
    Ts::Real = 0, minimal::Bool = false, contr::Bool = false, obs::Bool = false, noseig::Bool = false,
    atol::Real = zero(real(T)), rtol::Real = (min(size(P)...)*eps(real(float(one(T)))))*iszero(atol)) where T
    sysdata = pm2ls(P, contr = contr, obs = obs, noseig = noseig, minimal = minimal, atol = atol, rtol = rtol)
    return DescriptorStateSpace{eltype(sysdata[1])}(sysdata..., Ts)
end
dss(P::Union{VecOrMat{<:Polynomial},Polynomial}; kwargs...)  = dss(poly2pm(P); kwargs...)
function dss(P::LaurentPolynomial; kwargs...)
    n = degree(P)
    m = length(P)-1-n
    T = eltype(P)
    t = Polynomial{T,Polynomials.indeterminate(P)}([zeros(T,m);one(T)][:])
    dss(poly2pm(Polynomial(P*t)),poly2pm(t); kwargs...)
end
"""
    sys = dss(T, U, V, W; fast = true, contr = false, obs = false, minimal = false, atol = 0, rtol)

Construct an input-output equivalent descriptor system representation `sys = (A-λE,B,C,D)`
to a polynomial model specified by the polynomial matrices `T(λ)`, `U(λ)`, `V(λ)`, and `W(λ)`
such that

      V(λ)*inv(T(λ))*U(λ)+W(λ) = C*inv(λE-A)*B+D.

If `minimal = true`, the resulting realization `(A-λE,B,C,D)` has the least possible order `n` of `A-λE`.

`T(λ)`, `U(λ)`, `V(λ)`, and `W(λ)` can be specified as polynomial matrices of the form `X(λ) = X_1 + λ X_2 + ... + λ**k X_(k+1)`,
for `X = T`, `U`, `V`, and `W`, for which the coefficient matrices `X_i`, `i = 1, ..., k+1`, are stored in
the 3-dimensional matrices `X`, where `X[:,:,i]` contains the `i`-th coefficient matrix `X_i` (multiplying `λ**(i-1)`).

`T(λ)`, `U(λ)`, `V(λ)`, and `W(λ)` can also be specified as matrices, vectors or scalars of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `T(λ)`, `U(λ)`, `V(λ)`  and `W(λ)` have the same indeterminates.

The computed descriptor realization satisfies:

(1) `A-λE` is regular;

(2) `rank[B A-λE] = n` (controllability) if `minimal = true` or `contr = true`;

(3) `rank[A-λE; C] = n` (observability)  if `minimal = true` or `obs = true`;

(4) `A-λE` has no simple infinite eigenvalues if `minimal = true`.

The keyword arguments `atol` and `rtol`, specify, respectively, the absolute and relative tolerance for the
nonzero coefficients of the matrices `T(λ)`, `U(λ)`, `V(λ)` and `W(λ)`. The default relative tolerance is `nt*ϵ`,
where `nt` is the size of the square matrix `T(λ)` and `ϵ` is the machine epsilon of the element type of its coefficients.

The descriptor realization is built using the methods described in [1].

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).
"""
function dss(T::Array{T1,3}, U::Array{T2,3}, V::Array{T3,3}, W::Array{T4,3};
             Ts::Real = 0, contr::Bool = false, obs::Bool = false, minimal::Bool = false,
             fast::Bool = true, atol::Real = zero(real(T1)),
             rtol::Real = size(T,1)*eps(real(float(one(T1))))*iszero(atol)) where {T1, T2, T3, T4}
    sysdata = spm2ls(T, U, V, W, contr = contr, obs = obs, minimal = minimal, atol = atol, rtol = rtol)
    return DescriptorStateSpace{eltype(sysdata[1])}(sysdata..., Ts)
end
# function dss(T::Union{VecOrMat{Polynomial{T1}},VecOrMat{Polynomial{T1,X}},Polynomial{T1},Polynomial{T1,X}},
#              U::Union{VecOrMat{Polynomial{T2}},VecOrMat{Polynomial{T2,X}},Polynomial{T2},Polynomial{T2,X}},
#              V::Union{VecOrMat{Polynomial{T3}},VecOrMat{Polynomial{T3,X}},Polynomial{T3},Polynomial{T3,X}},
#              W::Union{VecOrMat{Polynomial{T4}},VecOrMat{Polynomial{T4,X}},Polynomial{T4},Polynomial{T4,X}}; kwargs...) where {T1, T2, T3, T4, X}
#              println("T = $T")
#        dss(poly2pm(T),poly2pm(U),poly2pm(V),poly2pm(W); kwargs...)
# end
# function dss(T::Union{VecOrMat{TP},TP},
#              U::Union{VecOrMat{TP},TP,Number,VecOrMat{<:Number}},
#              V::Union{VecOrMat{TP},TP,Number,VecOrMat{<:Number}},
#              W::Union{VecOrMat{TP},TP,Number,VecOrMat{<:Number}};
#              Ts::Real = 0, contr::Bool = false, obs::Bool = false, minimal::Bool = false,
#              fast::Bool = true, atol::Real = zero(float(real(eltype(eltype(T))))),
#              rtol::Real = size(T,1)*eps(one(eltype(atol)))*iszero(atol)) where {TP <: Polynomial}
#     return dss(poly2pm(T),poly2pm(U),poly2pm(V),poly2pm(W); Ts = Ts, contr = contr, obs = obs, minimal = minimal, fast = fast, atol = atol, rtol = rtol)
# end
function dss(T::Union{VecOrMat{<:Polynomial},Polynomial},
             U::Union{VecOrMat{<:Polynomial},Polynomial,Number,VecOrMat{<:Number}},
             V::Union{VecOrMat{<:Polynomial},Polynomial,Number,VecOrMat{<:Number}},
             W::Union{VecOrMat{<:Polynomial},Polynomial,Number,VecOrMat{<:Number}};
             Ts::Real = 0, contr::Bool = false, obs::Bool = false, minimal::Bool = false,
             fast::Bool = true, atol::Real = zero(float(real(eltype(eltype(T))))),
             rtol::Real = size(T,1)*eps(one(eltype(atol)))*iszero(atol))
return dss(poly2pm(T),poly2pm(U),poly2pm(V),poly2pm(W); Ts = Ts, contr = contr, obs = obs, minimal = minimal, fast = fast, atol = atol, rtol = rtol)
end

"""
    A, E, B, C, D  = dssdata([T,] sys)

Extract the matrices `A`, `E`, `B`, `C`, `D` of a descriptor system model `sys = (A-λE,B,C,D)`.
If the type `T` is specified, the resulting matrices are converted to this type.
"""
function dssdata(sys::DescriptorStateSpace)
    return copy(sys.A), copy(sys.E), copy(sys.B), copy(sys.C), copy(sys.D)
end
function dssdata(T::Type,sys::DescriptorStateSpace)
    #eltype(sys) == T ? (return sys.A, sys.E, sys.B, sys.C, sys.D) :
    return copy_oftype(sys.A,T), sys.E == I ? I : copy_oftype(sys.E,T), copy_oftype(sys.B,T), copy_oftype(sys.C,T), copy_oftype(sys.D,T)
end


text:
module WaterWaves1D

# external modules
using ProgressMeter,FFTW
using LinearMaps,IterativeSolvers,LinearAlgebra
using HDF5
import Base.show

# abstract types and structures
include("mesh.jl")
include("times.jl")
include("data.jl")
include("init.jl")
include("solver.jl")
include("models.jl")
include("problem.jl")

# tools
include("loadsave.jl")
include("tools.jl")
include("recipes.jl")

# initial data, models, solvers
include("initialdata/Random.jl")
include("initialdata/CnoidalWaveSerreGreenNaghdi.jl")
include("initialdata/SolitaryWaveSerreGreenNaghdi.jl")
include("initialdata/SolitaryWaveWhithamGreenNaghdi.jl")
include("initialdata/SolitaryWaveWhithamBoussinesq.jl")
include("initialdata/SolitaryWaveWhitham.jl")

include("models/Airy.jl")
include("models/Boussinesq.jl")
include("models/AkersNicholls.jl")
include("models/IsobeKakinuma.jl")
include("models/Matsuno.jl")
include("models/modifiedMatsuno.jl")
include("models/NonHydrostatic.jl")
include("models/SerreGreenNaghdi.jl")
include("models/SaintVenant.jl")
include("models/SquareRootDepth.jl")
include("models/WaterWaves.jl")
include("models/WhithamBoussinesq.jl")
include("models/WhithamGreenNaghdi.jl")
include("models/WWn.jl")

include("solvers/Euler.jl")
include("solvers/EulerSymp.jl")
include("solvers/RK4.jl")

end


text:
using GenTransportMCMC, Test

include("transport/transport.jl")

max_stars_repo_name: mvarble/gen-transport-mcmc

max_issues_repo_name: mvarble/gen-transport-mcmc

max_forks_repo_name: mvarble/gen-transport-mcmc


text:
using EEG, Lexicon

save("doc/API.md", EEG)


text:
module ImportAll

macro importall(moduleIn, verbose=false)
    modName = @eval($(moduleIn))
    out = []
    for name in names(modName)
            symbolname = @eval $(name)
            expression = :(import $(moduleIn): $(name))
            push!(out, expression)
    end

    ret = Expr(:block,out...)
	return ret
end

export @importall

end


text:
module Romania

export solveRomanianMapProblemBFS, solveRomanianMapProblemResultBFS,
       solveRomanianMapProblemUCS, solveRomanianMapProblemResultMinCost,
       solveRomanianMapProblemDLS, solveRomanianMapProblemResultDLS,
       solveRomanianMapProblemIDS, solveRomanianMapProblemResultIDS,
       solveRomanianMapProblemRBF, solveRomanianMapProblemResultRBF,
       solveRomanianMapProblemGSB, solveRomanianMapProblemResultGSB,
       solveRomanianMapProblemGSD, solveRomanianMapProblemResultGSD,
       solveRomanianMapProblemGSU, solveRomanianMapProblemResultGSU,
       solveRomanianMapProblemGSAS,solveRomanianMapProblemResultGSAS,
       solveRomanianMapProblemGSBF,solveRomanianMapProblemResultGSBF

using AIMACore

import Base: isempty, in
import AIMACore: search, goal_test, step_cost, actions, result, heuristic

struct In <: State
  place::Symbol
  In(place::Symbol) = new(place)
end
In(place::AbstractString) = In(Symbol(place))

struct Go <: Action
    place::Symbol
    Go(place::Symbol) = new(place)
end
Go(place::AbstractString) = Go(Symbol(place))

mutable struct RomaniaRoadMapProblem{SA <: SearchAlgorithm} <: Problem
    places::Dict{Symbol, Dict{Symbol,Int}}
    hSLD::Dict{Symbol, Int}
    initial_state::In
    goal_state::In
    search_algorithm::SA
    function RomaniaRoadMapProblem{SA}(init::In, goal::In,
                                       search_algorithm::SA) where {SA <: SearchAlgorithm}
        places = Dict(
            :Arad=>Dict(:Sibiu=>140,:Timisoara=>118,:Zerind=>75),
            :Bucharest=>Dict(:Fragaras=>211,:Giurgiu=>90,:Pitesti=>101,:Urziceni=>85),
            :Craiova=>Dict(:Drobeta=>120, :Pitesti=>138,:RimnicuVilcea=>146),
            :Drobeta=>Dict(:Craiova=>120, :Mehadia=>75),
            :Eforie=>Dict(:Hirsova=>86),
            :Fagaras=>Dict(:Bucharest=>211,:Sibiu=>99),
            :Giurgiu=>Dict(:Bucharest=>90),
            :Hirsova=>Dict(:Eforie=>86,:Urziceni=>98),
            :Iasi=>Dict(:Neamt=>87, :Vaslui=>92),
            :Lugoj=>Dict(:Mehadia=>70,:Timisoara=>111),
            :Mehadia=>Dict(:Drobeta=>75,:Lugoj=>70),
            :Neamt=>Dict(:Iasi=>87),
            :Oradea=>Dict(:Zerind=>71, :Sibiu=>151),
            :Pitesti=>Dict(:Bucharest=>101, :Craiova=>138, :RimnicuVilcea=>97),
            :RimnicuVilcea=>Dict(:Craiova=>146, :Pitesti=>97, :Sibiu=>80),
            :Sibiu=>Dict(:Arad=>140, :Fagaras=>99, :Oradea=>151, :RimnicuVilcea=>80),
            :Timisoara=>Dict(:Arad=>118, :Lugoj=>111),
            :Urziceni=>Dict(:Bucharest=>85, :Hirsova=>98, :Vaslui=>142),
            :Vaslui=>Dict(:Iasi=>92, :Urziceni=>142),
            :Zerind=>Dict(:Arad=>75, :Oradea=>71))

        hSLD = Dict(:Arad=>366,
                    :Bucharest=>0,
                    :Craiova=>160,
                    :Drobeta=>242,
                    :Eforie=>161,
                    :Fagaras=>176,
                    :Giurgiu=>77,
                    :Hirsova=>151,
                    :Iasi=>226,
                    :Lugoj=>244,
                    :Mehadia=>241,
                    :Neamt=>234,
                    :Oradea=>380,
                    :Pitesti=>100,
                    :RimnicuVilcea=>193,
                    :Sibiu=>253,
                    :Timisoara=>329,
                    :Urziceni=>80,
                    :Vaslui=>199,
                    :Zerind=>374)

        return new(places, hSLD, init, goal, search_algorithm)
    end
end

function result(problem::RomaniaRoadMapProblem, state::In, action::Go)
    s = state.place
    r = action.place
    haskey(problem.places[s], r) && return In(r)
    error("Invalid location cannot be reached")
end

function step_cost(problem::RomaniaRoadMapProblem, state::In, action::Go)
    s = state.place
    r = action.place
    return problem.places[s][r]
end

heuristic(problem::RomaniaRoadMapProblem, state::In) = problem.hSLD[state.place]

goal_test(problem::RomaniaRoadMapProblem, state::In) =
    (state == problem.goal_state)

function actions(problem::RomaniaRoadMapProblem, state::In)
    ret=[]
    for k in keys(problem.places[state.place])
        push!(ret, Go(k))
    end
    return ret
end

const solveRomanianMapProblemResultBFS =
    [(:Arad, 0), (:Sibiu, 140), (:Fagaras,239), (:Bucharest, 450)]
const solveRomanianMapProblemResultMinCost =
    [(:Arad, 0), (:Sibiu, 140), (:RimnicuVilcea, 220), (:Pitesti, 317),
     (:Bucharest, 418)]
const solveRomanianMapProblemResultDLS =
    [(:Arad, 0), (:Sibiu, 140), (:RimnicuVilcea, 220), (:Pitesti, 317),
     (:Bucharest, 418)]
const solveRomanianMapProblemResultIDS =
    [(:Arad, 0), (:Sibiu, 140), (:Fagaras,239), (:Bucharest, 450)]

const solveRomanianMapProblemResultGSB = solveRomanianMapProblemResultBFS
const solveRomanianMapProblemResultGSU = solveRomanianMapProblemResultMinCost
const solveRomanianMapProblemResultRBF = solveRomanianMapProblemResultMinCost
const solveRomanianMapProblemResultGSD = solveRomanianMapProblemResultIDS
const solveRomanianMapProblemResultGSBF = solveRomanianMapProblemResultIDS
const solveRomanianMapProblemResultGSAS = solveRomanianMapProblemResultMinCost

BFS = BreadthFirstSearch(In(""))
UCS = UniformCostSearch(In(""))
DLS = DepthLimitedSearch(4)
IDS = IterativeDeepeningSearch()
RBF = RecursiveBestFirstSearch(In(""))

GSB = GraphSearchBreadth(In(""))
GSD = GraphSearchDepth(In(""))
GSU = GraphSearchUniformCost(In(""))

GSBF = GraphSearchBestFirst(In(""))
GSAS = GraphSearchAStar(In(""))

function solveRomanianMapProblem(obj::T) where T
    problem = RomaniaRoadMapProblem{T}(In("Arad"), In("Bucharest"), obj)
    path = search(problem)
    path isa Symbol && return path
    ret=Tuple{Symbol, Int}[]
    for iter in path
        push!(ret, (iter.state.place, iter.path_cost))
    end
    return ret
end

solveRomanianMapProblemBFS() = solveRomanianMapProblem(BFS)
solveRomanianMapProblemUCS() = solveRomanianMapProblem(UCS)
solveRomanianMapProblemDLS() = solveRomanianMapProblem(DLS)
solveRomanianMapProblemIDS() = solveRomanianMapProblem(IDS)
solveRomanianMapProblemRBF() = solveRomanianMapProblem(RBF)
solveRomanianMapProblemGSB() = solveRomanianMapProblem(GSB)
solveRomanianMapProblemGSD() = solveRomanianMapProblem(GSD)
solveRomanianMapProblemGSU() = solveRomanianMapProblem(GSU)
solveRomanianMapProblemGSBF() = solveRomanianMapProblem(GSBF)
solveRomanianMapProblemGSAS() = solveRomanianMapProblem(GSAS)

end


text:
## perform the estimation of the mass segregation parameters in the e2e_products directory
## create the plots for each cluster
## output the results in a .csv file

rootdir= ENV["GAIA_ROOT"]
push!(LOAD_PATH,"$rootdir/run/src")
using GaiaClustering
using Printf, Glob, Statistics

import PyPlot , CSV , DataFrames

## directory
wdir    = "$rootdir/e2e_products"
plotdir = "$wdir/plotsIndividuals"
ocdir=    "$wdir/oc"

cd(wdir)

let
## Mass Segregation parameters#############################################
DTYPE= "2d"                 ## projection 2D/3D
DAVER= "geo"               ## average calculation: geo/ari
PERCMASSIVE= 20             ## percentile on the mag for massive/light stars

OUTFILE= "votlist-massSegregation.csv"
###########################################################################
@printf("## Compute the mass segregation parameters\n")
@printf("##\n")
############################functions #####################################
function plot_massSegregation(oc, figname, perc, kappa2d, kappa2derr)
    imass= select_massivestars(oc.gbar, perc, false)
    ilight= select_massivestars(oc.gbar, perc, true)

    xg= sum(oc.Y)/length(oc.Y) ;  yg= sum(oc.Z)/length(oc.Z)
    XX= oc.Y .- xg ; YY= oc.Z .- yg
    xmin= minimum(XX) ; xmax= maximum(XX)
    ymin= minimum(YY) ; ymax= maximum(YY)
    dl= max(xmax-xmin,ymax-ymin)*1.05 ; x2= (xmin+xmax)/2 ; y2= (ymin+ymax)/2

    fig= PyPlot.figure(figsize=(7,7))
    ax= PyPlot.subplot(111)
    PyPlot.grid("on")

    rr= sqrt.((XX .* XX) .+ (YY .* YY))
    siglight= std(rr[ilight]) ; sigmass= std(rr[imass])

    circle = PyPlot.plt.Circle((0,0), siglight, color="b", fill=false, linestyle="--")
    ax.add_artist(circle)
    circle = PyPlot.plt.Circle((0,0), sigmass, color="r", fill=false, linestyle="--")
    ax.add_artist(circle)

    PyPlot.title(figname[1:end-20])
    PyPlot.xlabel("Y(pc)") ; PyPlot.ylabel("Z(pc)")
    PyPlot.xlim([x2-dl/2,x2+dl/2]) ; PyPlot.ylim([y2-dl/2,y2+dl/2])

    PyPlot.plt.plot(XX, YY,".k", ms= 3.0, label="all")
    PyPlot.plt.plot(XX[imass], YY[imass],"or",label="massive")
    PyPlot.plt.plot(XX[ilight], YY[ilight],"ob",label="light")
    PyPlot.legend()
    txt= @sprintf("kappa: %3.2f(%3.2f)",kappa2d,kappa2derr)
    PyPlot.text(x2-dl*0.45,y2-dl*0.45,txt)
    PyPlot.savefig(figname)
end
## main loop
    cd(ocdir)
    files= glob("*-oc.csv")
    cd(wdir)

    ilabel= []
    Κms2d= Array{Float64,1}(undef,0) ; Κms3d = Array{Float64,1}(undef,0) ; k3=  Array{Float64,1}(undef,0)
    Q=  Array{Float64,1}(undef,0)
    errΚms2d= Array{Float64,1}(undef,0) ; errΚms3d= Array{Float64,1}(undef,0) ; errk3= Array{Float64,1}(undef,0)
    votname= Array{String,1}(undef,0)

    i= 0
    for f in files
        i+= 1
        oc= CSV.read("$ocdir/$f" , delim= ";")

        vot= @sprintf("%s.vot", f[1:end-7])
        println(vot)
        push!(votname, vot )

        y, yerr= get_kappaMS(oc[:gbar],oc[:X], oc[:Y] , oc[:Z],PERCMASSIVE , "3d", DAVER, false)
        push!(ilabel, i)
        push!(Κms3d, y)
        push!(errΚms3d, yerr)

        y, yerr= get_kappaMS(oc[:gbar],oc[:X], oc[:Y] , oc[:Z], PERCMASSIVE, "2d", DAVER, false)
        push!(Κms2d, y)
        push!(errΚms2d, yerr)

        y, yerr= get_kappaMS(oc[:gbar],oc[:vl], oc[:vl] , oc[:vb], PERCMASSIVE, "2d", DAVER, false)
        push!(k3, y)
        push!(errk3, yerr)

        y= get_Q(oc[:X], oc[:Y] , oc[:Z])
        push!(Q, y)

        cd(plotdir)
        figname= @sprintf("%s-massSegregation.png", f[1:end-7])
        plot_massSegregation(oc, figname, PERCMASSIVE, Κms2d[end],errΚms2d[end])
        cd(wdir)
    end

    df= DataFrames.DataFrame()
    df.votname= votname
    df.kappa2d= Κms2d ; df.kappa2derr= errΚms2d
    df.kappa3d= Κms3d ; df.kappa3derr= errΚms3d
    df.kappavel= k3 ; df.kappavelerr= errk3
    df.Q= Q

    CSV.write(OUTFILE, df, delim=";")
    println(df)

    @printf("### kappa2d median: %3.3f(%3.3f)\n",median(Κms2d),std(Κms2d))
    @printf("### kappa3d median: %3.3f(%3.3f)\n",median(Κms3d),std(Κms3d))
    @printf("### kappa_vel median: %3.3f(%3.3f)\n",median(k3),std(k3))
    @printf("### Q(clustering) median: %3.3f(%3.3f)\n",median(Q),std(Q))
end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
using AMDGPUnative
const AS = AMDGPUnative.AS
using GPUCompiler
using LLVM, LLVM.Interop
using HSARuntime
using InteractiveUtils
using SpecialFunctions
using Test

agent_name = HSARuntime.get_name(get_default_agent())
agent_isa = get_first_isa(get_default_agent())
@info "Testing using device $agent_name with ISA $agent_isa"

@testset "AMDGPUnative" begin

@testset "Core" begin
include("pointer.jl")
end

if AMDGPUnative.configured
    @test length(get_agents()) > 0
    if length(get_agents()) > 0
        include("synchronization.jl")
        include("trap.jl")
        @testset "Device" begin
            include("device/launch.jl")
            include("device/vadd.jl")
            include("device/memory.jl")
            include("device/indexing.jl")
            include("device/hostcall.jl")
            include("device/output.jl")
            include("device/globals.jl")
            if Base.libllvm_version >= v"7.0"
                include("device/math.jl")
            else
                @warn "Testing with LLVM 6; some tests will be disabled!"
                @test_skip "Math Intrinsics"
            end
        end
    end
else
    @warn("AMDGPUnative.jl has not been configured; skipping on-device tests.")
end

end


text:
import UUIDs

# This function is taken from:
# https://github.com/JuliaRegistries/RegistryCI.jl/blob/master/src/RegistryCI.jl
function gather_stdlib_uuids()
    if VERSION < v"1.1"
        return Set{UUIDs.UUID}(x for x in keys(Pkg.Types.gather_stdlib_uuids()))
    else
        return Set{UUIDs.UUID}(x for x in keys(Pkg.Types.stdlib()))
    end
end


text:
import Pkg; Pkg.add("DataFramesMeta")
using Random
using LinearAlgebra
using Statistics
using Optim
using DataFrames
using DataFramesMeta
using CSV
using HTTP
using GLM
# read in function to create state transitions for dynamic model
#The function didn't work for me for some reason - I tried it a couple of times but couldn't resolve the error so I'm copying the code from the file you uploaded which creates the grid#
function create_grids()

    function xgrid(theta,xval)
        N      = length(xval)
        xub    = vcat(xval[2:N],Inf)
        xtran1  = zeros(N,N)
        xtran1c = zeros(N,N)
        lcdf   = zeros(N)
        for i=1:length(xval)
            xtran1[:,i]   = (xub[i] .>= xval).*(1 .- exp.(-theta*(xub[i] .- xval)) .- lcdf)
            lcdf        .+= xtran1[:,i]
            xtran1c[:,i] .+= lcdf
        end
        return xtran1, xtran1c
    end

    zval   = collect([.25:.01:1.25]...)
    zbin   = length(zval)
    xval   = collect([0:.125:25]...)
    xbin   = length(xval)
    tbin   = xbin*zbin
    xtran  = zeros(tbin,xbin)
    xtranc = zeros(xbin,xbin,xbin)
    for z=1:zbin
        xtran[(z-1)*xbin+1:z*xbin,:],xtranc[:,:,z] = xgrid(zval[z],xval)
    end

    return zval,zbin,xval,xbin,xtran
end

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 1: reshaping the data
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# load in the data
using DataFrames
using CSV
using HTTP
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS5-ddc/busdataBeta0.csv"
df = CSV.read(HTTP.get(url).body, DataFrame)
# create bus id variable
df = @transform(df, bus_id = 1:size(df,1))

#---------------------------------------------------
# reshape from wide to long (must do this twice be-
# cause DataFrames.stack() requires doing it one
# variable at a time)
#---------------------------------------------------
# first reshape the decision variable
dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

# next reshape the odometer variable
dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

# join reshaped df's back together
df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 2: estimate a static version of the model using GLM
#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
theta_hat = glm(@formula(Y ~ Odometer + Branded ), df_long, Binomial(), LogitLink())
println(theta_hat)
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3a,b,c: read in data for dynamic model
#:::::::::::::::::::::::::::::::::::::::::::::::::::
using DataFrames
using CSV
using HTTP
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS5-ddc/busdata.csv"
df = CSV.read(HTTP.get(url).body, DataFrame)
Y = Matrix(df[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])
# [ convert other data frame columns to matrices]
zval,zbin,xval,xbin,xtran = create_grids()
@views @inbounds function myfun(theta)   #3f?#
N=size(xtran,1) #
T=size(Y,2) #
FV=zeros(N,2,T+1) #
#Now write four nested for loops over each of the possible states#
    for t in T:-1:1 #Loop backwards over t from T +1 to 1#
    for b in 0:1#Loop over the two possible brand states {0,1}#
    for z in 1:zbin#Loop over the possible permanent route usage states (i.e. from 1 to zbin)#
    for x in 1:xbin#Loop over the possible odometer states (i.e. from 1 to xbin#
    rowtm=x+(z-1)*xbin #Create an object that marks the row of the transition matrix that we need to - This will be x + (z-1)*xbin (where x indexes the mileage bin and z indexes the route usage bin), given how the xtran matrix was constructed in the create grids() function.#
    end
    #conditional value function#
v1=theta[1]+theta[2]*xval[x]+theta[3]*b+ beta*xtran[row,:]'*FV[(z-1)*xbin+1:z*xbin,b+1,t+1]
v0=beta*xtran[1+(z-1)*xbin,:]'*FV[(z-1)*xbin+1:z*xbin,b+1,t+1]
FV[row,b+1,t]=beta*log(exp(v1)+exp(v0))
beta=0.9
#3d#
loglikelihood=0
v_0=beta*xtran[1+(:Zst-1)*xbin,:]'*FV[(:Zst-1)*xbin+1:Zst*xbin,b+1,t+1]
v_1=theta[1]+theta[2]*xval[:xst]+theta[3]*b+ beta*xtran[row,:]'*FV[(:zst-1)*xbin+1:z*xbin,b+1,t+1]
#3e#
return=-loglikelihood


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

struct LearningSolver
    py::PyCall.PyObject
    optimizer_factory::Any
end


text:
#
# @description Unit tests for PKCS#7 padding method.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2015 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
@testset "PKCS#7 padder tests" begin
  padder = Pkcs7Padder()

  @test pad(padder, [ 0x0A, 0x0B ], 8) == [ 0x0A, 0x0B, 0x06, 0x06, 0x06, 0x06, 0x06, 0x06 ]
  @test pad(padder, [ 0x0A, 0x0B ], 2) == [ 0x0A, 0x0B, 0x02, 0x02 ]

  @test unpad(padder, [ 0x0A, 0x0B, 0x06, 0x06, 0x06, 0x06, 0x06, 0x06 ]) == [ 0x0A, 0x0B ]
  @test unpad(padder, [ 0x02, 0x02 ]) == []

  @test_throws PaddingError unpad(padder, [ 0xFF, 0x00 ])
  @test_throws PaddingError unpad(padder, [ 0x03, 0x03 ])
  @test_throws PaddingError unpad(padder, [ 0x0A, 0x0B, 0x03, 0x03 ])
  @test_throws PaddingError unpad(padder, [ 0x0A, 0x0B, 0x0F, 0x02 ])
end


text:
const  √ = sqrt
const  ∑ = sum
const  ⊗ = kron

function hom⁻¹(v::StaticVector)
    if isapprox(v[end], 0.0; atol = 1e-14)
        pop(v)
    else
        pop(v / v[end])
    end
end

function hom(v::StaticVector)
    push(v,1)
end

# Used to convert a N x 2 matrix into a one-dimensional array of 2D points (each point is a row of the original matrix).
function svectors(x::AbstractArray{T}, ::Val{N}) where {T,N}
    size(x,1) == N || error("sizes mismatch")
    isbitstype(T) || error("use for bitstypes only")
    reinterpret(SVector{N,T}, vec(x))
end


text:
println("Using the Task's Test Data")
test = """12 127 28 42 39 113 42 18 44 118 44 37 113 124 37 48 127 36 29
   31 125 139 131 115 105 132 104 123 35 113 122 42 117 119 58 109 23 105
   63 27 44 105 99 41 128 121 116 125 32 61 37 127 29 113 121 58 114 126
   53 114 96 25 109 7 31 141 46 13 27 43 117 116 27 7 68 40 31 115 124 42
   128 52 71 118 117 38 27 106 33 117 116 111 40 119 47 105 57 122 109
   124 115 43 120 43 27 27 18 28 48 125 107 114 34 133 45 120 30 127 31
   116 146"""
test = map(parseint, split(test, r"\s"))
println(stemleaf(test))

println("Test with Reals and Negative Zero Stem")
test = [-23.678758, -12.45, -3.4, 4.43, 5.5, 5.678, 16.87, 24.7, 56.8]
println(stemleaf(test))

println("Test with Leaf Size Scaling")
test = int(500*randn(20))
println("Using:  ", test)
println(stemleaf(test, 10))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
Implement a queue using two stacks. Recall that a queue is a FIFO (first-in, first-out) data structure with the following methods: enqueue, which inserts an element into the queue, and dequeue, which removes it.
=#

mutable struct Queue
    stack1::Array{Any, 1}
    stack2::Array{Any, 1}
end

function enqueue(queue::Queue, element::Any)
    push!(queue.stack1, element)
end

function dequeue(queue::Queue)
    while length(queue.stack1) > 0
        push!(queue.stack2, pop!(queue.stack1))
    end
    out = pop!(queue.stack2)
    while length(queue.stack2) > 0
        push!(queue.stack1, pop!(queue.stack2))
    end
    return out
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem