text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                                                                  %%%%%
%%%%    IEEE PES Power Grid Library - Optimal Power Flow - v21.07     %%%%%
%%%%          (https://github.com/power-grid-lib/pglib-opf)           %%%%%
%%%%             Benchmark Group - Active Power Increase              %%%%%
%%%%                         29 - July - 2021                         %%%%%
%%%%                                                                  %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mpc = pglib_opf_case39_epri__api
mpc.version = '2';
mpc.baseMVA = 100.0;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [
	1	 1	 160.91	 44.20	 0.0	 0.0	 2	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	2	 1	 0.0	 0.0	 0.0	 0.0	 2	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	3	 1	 530.88	 2.40	 0.0	 0.0	 2	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	4	 1	 824.35	 184.00	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	5	 1	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	6	 1	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	7	 1	 385.47	 84.00	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	8	 1	 860.62	 176.60	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	9	 1	 6.50	 -66.60	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	10	 1	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	11	 1	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	12	 1	 14.06	 88.00	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	13	 1	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	14	 1	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	15	 1	 527.58	 153.00	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	16	 1	 542.42	 32.30	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	17	 1	 0.0	 0.0	 0.0	 0.0	 2	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	18	 1	 260.49	 30.00	 0.0	 0.0	 2	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	19	 1	 0.0	 0.0	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	20	 1	 1121.12	 103.00	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	21	 1	 451.74	 115.00	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	22	 1	 0.0	 0.0	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	23	 1	 408.05	 84.60	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	24	 1	 308.60	 -92.20	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	25	 1	 369.31	 47.20	 0.0	 0.0	 2	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	26	 1	 229.17	 17.00	 0.0	 0.0	 2	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	27	 1	 463.28	 75.50	 0.0	 0.0	 2	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	28	 1	 339.63	 27.60	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	29	 1	 467.41	 26.90	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	30	 2	 0.0	 0.0	 0.0	 0.0	 2	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	31	 3	 15.17	 4.60	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	32	 2	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	33	 2	 0.0	 0.0	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	34	 2	 0.0	 0.0	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	35	 2	 0.0	 0.0	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	36	 2	 0.0	 0.0	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	37	 2	 0.0	 0.0	 0.0	 0.0	 2	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	38	 2	 0.0	 0.0	 0.0	 0.0	 3	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
	39	 2	 1820.16	 250.00	 0.0	 0.0	 1	    1.00000	    0.00000	 345.0	 1	    1.06000	    0.94000;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin
mpc.gen = [
	30	 159.0	 59.0	 400.0	 -282.0	 1.0	 100.0	 1	 318	 0.0; % COW
	31	 789.0	 0.0	 854.4	 -854.4	 1.0	 100.0	 1	 1578	 0.0; % COW
	32	 1277.5	 0.0	 1278.0	 -1278.0	 1.0	 100.0	 1	 2555	 0.0; % COW
	33	 432.5	 0.0	 433.0	 -433.0	 1.0	 100.0	 1	 865	 0.0; % COW
	34	 438.5	 0.0	 439.0	 -439.0	 1.0	 100.0	 1	 877	 0.0; % COW
	35	 560.0	 0.0	 560.0	 -560.0	 1.0	 100.0	 1	 1120	 0.0; % COW
	36	 461.0	 0.0	 461.0	 -461.0	 1.0	 100.0	 1	 922	 0.0; % COW
	37	 759.0	 0.0	 759.0	 -759.0	 1.0	 100.0	 1	 1518	 0.0; % COW
	38	 594.5	 0.0	 595.0	 -595.0	 1.0	 100.0	 1	 1189	 0.0; % COW
	39	 1261.5	 0.0	 1262.0	 -1262.0	 1.0	 100.0	 1	 2523	 0.0; % COW
];

%% generator cost data
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
	2	 0.0	 0.0	 3	   0.000000	   6.724778	   0.000000; % COW
	2	 0.0	 0.0	 3	   0.000000	  14.707625	   0.000000; % COW
	2	 0.0	 0.0	 3	   0.000000	  24.804734	   0.000000; % COW
	2	 0.0	 0.0	 3	   0.000000	  34.844643	   0.000000; % COW
	2	 0.0	 0.0	 3	   0.000000	  24.652994	   0.000000; % COW
	2	 0.0	 0.0	 3	   0.000000	  32.306483	   0.000000; % COW
	2	 0.0	 0.0	 3	   0.000000	  18.157477	   0.000000; % COW
	2	 0.0	 0.0	 3	   0.000000	  31.550181	   0.000000; % COW
	2	 0.0	 0.0	 3	   0.000000	  22.503168	   0.000000; % COW
	2	 0.0	 0.0	 3	   0.000000	  27.434444	   0.000000; % COW
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [
	1	 2	 0.0035	 0.0411	 0.6987	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	1	 39	 0.001	 0.025	 0.75	 1000.0	 1000.0	 1000.0	 0.0	 0.0	 1	 -30.0	 30.0;
	2	 3	 0.0013	 0.0151	 0.2572	 500.0	 500.0	 500.0	 0.0	 0.0	 1	 -30.0	 30.0;
	2	 25	 0.007	 0.0086	 0.146	 500.0	 500.0	 500.0	 0.0	 0.0	 1	 -30.0	 30.0;
	2	 30	 0.0	 0.0181	 0.0	 900.0	 900.0	 2500.0	 1.025	 0.0	 1	 -30.0	 30.0;
	3	 4	 0.0013	 0.0213	 0.2214	 500.0	 500.0	 500.0	 0.0	 0.0	 1	 -30.0	 30.0;
	3	 18	 0.0011	 0.0133	 0.2138	 500.0	 500.0	 500.0	 0.0	 0.0	 1	 -30.0	 30.0;
	4	 5	 0.0008	 0.0128	 0.1342	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	4	 14	 0.0008	 0.0129	 0.1382	 500.0	 500.0	 500.0	 0.0	 0.0	 1	 -30.0	 30.0;
	5	 6	 0.0002	 0.0026	 0.0434	 1200.0	 1200.0	 1200.0	 0.0	 0.0	 1	 -30.0	 30.0;
	5	 8	 0.0008	 0.0112	 0.1476	 900.0	 900.0	 900.0	 0.0	 0.0	 1	 -30.0	 30.0;
	6	 7	 0.0006	 0.0092	 0.113	 900.0	 900.0	 900.0	 0.0	 0.0	 1	 -30.0	 30.0;
	6	 11	 0.0007	 0.0082	 0.1389	 480.0	 480.0	 480.0	 0.0	 0.0	 1	 -30.0	 30.0;
	6	 31	 0.0	 0.025	 0.0	 1800.0	 1800.0	 1800.0	 1.07	 0.0	 1	 -30.0	 30.0;
	7	 8	 0.0004	 0.0046	 0.078	 900.0	 900.0	 900.0	 0.0	 0.0	 1	 -30.0	 30.0;
	8	 9	 0.0023	 0.0363	 0.3804	 900.0	 900.0	 900.0	 0.0	 0.0	 1	 -30.0	 30.0;
	9	 39	 0.001	 0.025	 1.2	 900.0	 900.0	 900.0	 0.0	 0.0	 1	 -30.0	 30.0;
	10	 11	 0.0004	 0.0043	 0.0729	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	10	 13	 0.0004	 0.0043	 0.0729	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	10	 32	 0.0	 0.02	 0.0	 900.0	 900.0	 2500.0	 1.07	 0.0	 1	 -30.0	 30.0;
	12	 11	 0.0016	 0.0435	 0.0	 500.0	 500.0	 500.0	 1.006	 0.0	 1	 -30.0	 30.0;
	12	 13	 0.0016	 0.0435	 0.0	 500.0	 500.0	 500.0	 1.006	 0.0	 1	 -30.0	 30.0;
	13	 14	 0.0009	 0.0101	 0.1723	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	14	 15	 0.0018	 0.0217	 0.366	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	15	 16	 0.0009	 0.0094	 0.171	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	16	 17	 0.0007	 0.0089	 0.1342	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	16	 19	 0.0016	 0.0195	 0.304	 600.0	 600.0	 2500.0	 0.0	 0.0	 1	 -30.0	 30.0;
	16	 21	 0.0008	 0.0135	 0.2548	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	16	 24	 0.0003	 0.0059	 0.068	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	17	 18	 0.0007	 0.0082	 0.1319	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	17	 27	 0.0013	 0.0173	 0.3216	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	19	 20	 0.0007	 0.0138	 0.0	 900.0	 900.0	 2500.0	 1.06	 0.0	 1	 -30.0	 30.0;
	19	 33	 0.0007	 0.0142	 0.0	 900.0	 900.0	 2500.0	 1.07	 0.0	 1	 -30.0	 30.0;
	20	 34	 0.0009	 0.018	 0.0	 900.0	 900.0	 2500.0	 1.009	 0.0	 1	 -30.0	 30.0;
	21	 22	 0.0008	 0.014	 0.2565	 900.0	 900.0	 900.0	 0.0	 0.0	 1	 -30.0	 30.0;
	22	 23	 0.0006	 0.0096	 0.1846	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	22	 35	 0.0	 0.0143	 0.0	 900.0	 900.0	 2500.0	 1.025	 0.0	 1	 -30.0	 30.0;
	23	 24	 0.0022	 0.035	 0.361	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	23	 36	 0.0005	 0.0272	 0.0	 900.0	 900.0	 2500.0	 0.0	 0.0	 1	 -30.0	 30.0;
	25	 26	 0.0032	 0.0323	 0.531	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	25	 37	 0.0006	 0.0232	 0.0	 900.0	 900.0	 2500.0	 1.025	 0.0	 1	 -30.0	 30.0;
	26	 27	 0.0014	 0.0147	 0.2396	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	26	 28	 0.0043	 0.0474	 0.7802	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	26	 29	 0.0057	 0.0625	 1.029	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	28	 29	 0.0014	 0.0151	 0.249	 600.0	 600.0	 600.0	 0.0	 0.0	 1	 -30.0	 30.0;
	29	 38	 0.0008	 0.0156	 0.0	 1200.0	 1200.0	 2500.0	 1.025	 0.0	 1	 -30.0	 30.0;
];

% INFO    : === Translation Options ===
% INFO    : Load Model:                  from file ./pglib_opf_case39_epri.m.api.sol
% INFO    : Gen Active Capacity Model:   stat
% INFO    : Gen Reactive Capacity Model: al50ag
% INFO    : Gen Active Cost Model:       stat
% INFO    :
% INFO    : === Load Replacement Notes ===
% INFO    : Bus 1	: Pd=97.6, Qd=44.2 -> Pd=160.91, Qd=44.20
% INFO    : Bus 3	: Pd=322.0, Qd=2.4 -> Pd=530.88, Qd=2.40
% INFO    : Bus 4	: Pd=500.0, Qd=184.0 -> Pd=824.35, Qd=184.00
% INFO    : Bus 7	: Pd=233.8, Qd=84.0 -> Pd=385.47, Qd=84.00
% INFO    : Bus 8	: Pd=522.0, Qd=176.6 -> Pd=860.62, Qd=176.60
% INFO    : Bus 9	: Pd=6.5, Qd=-66.6 -> Pd=6.50, Qd=-66.60
% INFO    : Bus 12	: Pd=8.53, Qd=88.0 -> Pd=14.06, Qd=88.00
% INFO    : Bus 15	: Pd=320.0, Qd=153.0 -> Pd=527.58, Qd=153.00
% INFO    : Bus 16	: Pd=329.0, Qd=32.3 -> Pd=542.42, Qd=32.30
% INFO    : Bus 18	: Pd=158.0, Qd=30.0 -> Pd=260.49, Qd=30.00
% INFO    : Bus 20	: Pd=680.0, Qd=103.0 -> Pd=1121.12, Qd=103.00
% INFO    : Bus 21	: Pd=274.0, Qd=115.0 -> Pd=451.74, Qd=115.00
% INFO    : Bus 23	: Pd=247.5, Qd=84.6 -> Pd=408.05, Qd=84.60
% INFO    : Bus 24	: Pd=308.6, Qd=-92.2 -> Pd=308.60, Qd=-92.20
% INFO    : Bus 25	: Pd=224.0, Qd=47.2 -> Pd=369.31, Qd=47.20
% INFO    : Bus 26	: Pd=139.0, Qd=17.0 -> Pd=229.17, Qd=17.00
% INFO    : Bus 27	: Pd=281.0, Qd=75.5 -> Pd=463.28, Qd=75.50
% INFO    : Bus 28	: Pd=206.0, Qd=27.6 -> Pd=339.63, Qd=27.60
% INFO    : Bus 29	: Pd=283.5, Qd=26.9 -> Pd=467.41, Qd=26.90
% INFO    : Bus 31	: Pd=9.2, Qd=4.6 -> Pd=15.17, Qd=4.60
% INFO    : Bus 39	: Pd=1104.0, Qd=250.0 -> Pd=1820.16, Qd=250.00
% INFO    :
% INFO    : === Generator Setpoint Replacement Notes ===
% INFO    : Gen at bus 30	: Pg=520.0, Qg=270.0 -> Pg=90.0, Qg=235.0
% INFO    : Gen at bus 31	: Pg=323.0, Qg=100.0 -> Pg=1455.0, Qg=712.0
% INFO    : Gen at bus 32	: Pg=362.5, Qg=225.0 -> Pg=805.0, Qg=402.0
% INFO    : Gen at bus 33	: Pg=326.0, Qg=125.0 -> Pg=858.0, Qg=253.0
% INFO    : Gen at bus 34	: Pg=254.0, Qg=83.5 -> Pg=859.0, Qg=251.0
% INFO    : Gen at bus 35	: Pg=343.5, Qg=100.0 -> Pg=864.0, Qg=250.0
% INFO    : Gen at bus 36	: Pg=290.0, Qg=120.0 -> Pg=874.0, Qg=215.0
% INFO    : Gen at bus 37	: Pg=282.0, Qg=125.0 -> Pg=887.0, Qg=3.0
% INFO    : Gen at bus 38	: Pg=432.5, Qg=75.0 -> Pg=1185.0, Qg=188.0
% INFO    : Gen at bus 39	: Pg=550.0, Qg=100.0 -> Pg=2326.0, Qg=143.0
% INFO    :
% INFO    : === Generator Reactive Capacity Atleast Setpoint Value Notes ===
% INFO    : Gen at bus 30	: Qg 235.0, Qmin 140.0, Qmax 400.0 -> Qmin -282.0, Qmax 400.0
% INFO    : Gen at bus 31	: Qg 712.0, Qmin -100.0, Qmax 300.0 -> Qmin -854.4, Qmax 854.4
% INFO    : Gen at bus 32	: Qg 402.0, Qmin 150.0, Qmax 300.0 -> Qmin -482.4, Qmax 482.4
% INFO    : Gen at bus 33	: Qg 253.0, Qmin 0.0, Qmax 250.0 -> Qmin -303.6, Qmax 303.6
% INFO    : Gen at bus 34	: Qg 251.0, Qmin 0.0, Qmax 167.0 -> Qmin -301.2, Qmax 301.2
% INFO    : Gen at bus 35	: Qg 250.0, Qmin -100.0, Qmax 300.0 -> Qmin -300.0, Qmax 300.0
% INFO    : Gen at bus 36	: Qg 215.0, Qmin 0.0, Qmax 240.0 -> Qmin -258.0, Qmax 240.0
% INFO    : Gen at bus 37	: Qg 3.0, Qmin 0.0, Qmax 250.0 -> Qmin -3.6, Qmax 250.0
% INFO    : Gen at bus 38	: Qg 188.0, Qmin -150.0, Qmax 300.0 -> Qmin -225.6, Qmax 300.0
% INFO    : Gen at bus 39	: Qg 143.0, Qmin -100.0, Qmax 300.0 -> Qmin -171.6, Qmax 300.0
% INFO    :
% INFO    : === Generator Classification Notes ===
% INFO    : COW    10  -   100.00
% INFO    :
% INFO    : === Generator Active Capacity Stat Model Notes ===
% INFO    : Gen at bus 30 - COW	: Pg=90.0, Pmax=1040.0 -> Pmax=318   samples: 2
% WARNING : Failed to find a generator capacity within (1455.0-7275.0) after 100 samples, using percent increase model
% INFO    : Gen at bus 31 - COW	: Pg=1455.0, Pmax=646.0 -> Pmax=1578   samples: 100
% INFO    : Gen at bus 32 - COW	: Pg=805.0, Pmax=725.0 -> Pmax=2555   samples: 26
% INFO    : Gen at bus 33 - COW	: Pg=858.0, Pmax=652.0 -> Pmax=865   samples: 44
% INFO    : Gen at bus 34 - COW	: Pg=859.0, Pmax=508.0 -> Pmax=877   samples: 33
% INFO    : Gen at bus 35 - COW	: Pg=864.0, Pmax=687.0 -> Pmax=1120   samples: 15
% INFO    : Gen at bus 36 - COW	: Pg=874.0, Pmax=580.0 -> Pmax=922   samples: 12
% INFO    : Gen at bus 37 - COW	: Pg=887.0, Pmax=564.0 -> Pmax=1518   samples: 39
% INFO    : Gen at bus 38 - COW	: Pg=1185.0, Pmax=865.0 -> Pmax=1189   samples: 57
% WARNING : Failed to find a generator capacity within (2326.0-11630.0) after 100 samples, using percent increase model
% INFO    : Gen at bus 39 - COW	: Pg=2326.0, Pmax=1100.0 -> Pmax=2523   samples: 100
% INFO    :
% INFO    : === Generator Active Capacity LB Model Notes ===
% INFO    :
% INFO    : === Generator Reactive Capacity Atleast Max 50 Percent Active Model Notes ===
% INFO    : Gen at bus 32 - COW	: Pmax 2555.0, Qmin -482.4, Qmax 482.4 -> Qmin -1278.0, Qmax 1278.0
% INFO    : Gen at bus 33 - COW	: Pmax 865.0, Qmin -303.6, Qmax 303.6 -> Qmin -433.0, Qmax 433.0
% INFO    : Gen at bus 34 - COW	: Pmax 877.0, Qmin -301.2, Qmax 301.2 -> Qmin -439.0, Qmax 439.0
% INFO    : Gen at bus 35 - COW	: Pmax 1120.0, Qmin -300.0, Qmax 300.0 -> Qmin -560.0, Qmax 560.0
% INFO    : Gen at bus 36 - COW	: Pmax 922.0, Qmin -258.0, Qmax 240.0 -> Qmin -461.0, Qmax 461.0
% INFO    : Gen at bus 37 - COW	: Pmax 1518.0, Qmin -3.6, Qmax 250.0 -> Qmin -759.0, Qmax 759.0
% INFO    : Gen at bus 38 - COW	: Pmax 1189.0, Qmin -225.6, Qmax 300.0 -> Qmin -595.0, Qmax 595.0
% INFO    : Gen at bus 39 - COW	: Pmax 2523.0, Qmin -171.6, Qmax 300.0 -> Qmin -1262.0, Qmax 1262.0
% INFO    :
% INFO    : === Generator Setpoint Replacement Notes ===
% INFO    : Gen at bus 30	: Pg=90.0, Qg=235.0 -> Pg=159.0, Qg=59.0
% INFO    : Gen at bus 30	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 31	: Pg=1455.0, Qg=712.0 -> Pg=789.0, Qg=0.0
% INFO    : Gen at bus 31	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 32	: Pg=805.0, Qg=402.0 -> Pg=1277.5, Qg=0.0
% INFO    : Gen at bus 32	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 33	: Pg=858.0, Qg=253.0 -> Pg=432.5, Qg=0.0
% INFO    : Gen at bus 33	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 34	: Pg=859.0, Qg=251.0 -> Pg=438.5, Qg=0.0
% INFO    : Gen at bus 34	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 35	: Pg=864.0, Qg=250.0 -> Pg=560.0, Qg=0.0
% INFO    : Gen at bus 35	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 36	: Pg=874.0, Qg=215.0 -> Pg=461.0, Qg=0.0
% INFO    : Gen at bus 36	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 37	: Pg=887.0, Qg=3.0 -> Pg=759.0, Qg=0.0
% INFO    : Gen at bus 37	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 38	: Pg=1185.0, Qg=188.0 -> Pg=594.5, Qg=0.0
% INFO    : Gen at bus 38	: Vg=1.0 -> Vg=1.0
% INFO    : Gen at bus 39	: Pg=2326.0, Qg=143.0 -> Pg=1261.5, Qg=0.0
% INFO    : Gen at bus 39	: Vg=1.0 -> Vg=1.0
% INFO    :
% INFO    : === Writing Matpower Case File Notes ===

author: power-grid-lib

repo: pglib-opf


text:
classdef FuncHelpClass < matlab.mixin.Copyable

    properties
        funcname;
        helpstr;
        sections;
        userfuncdir;
        config;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  Initialization methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % -------------------------------------------------------------
        function obj = FuncHelpClass(funcname)
            % Syntax:
            %       fhelp = FuncHelpClass(funcname);
            %
            % Examples:
            %       fhelp = FuncHelpClass('hmrR_DeconvHRF_DriftSS');
            %
            obj.userfuncdir = {};
            obj.config = [];

            if nargin==0
                return;
            end
            obj.funcname = funcname;
            if isdeployed()
                obj.userfuncdir = obj.UserFuncDirs(obj);
            end
            obj.Help();
            obj.ParseSections();

        end

        % ---------------------------------------------------------------
        function AddSection(obj, name)
            fieldname = name;
            fieldname(fieldname==' ')='';
            fieldname = lower(fieldname);
            newsection = struct(...
                'name',name, ...
                'lines',[0,0], ...
                'str','', ...
                'subsections', [] ...
                );
            eval( sprintf('obj.sections.%s = newsection;', fieldname) );
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  Methods for parsing sections
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % -------------------------------------------------------------
        function ParseSections(obj)
            % This method dynamically parses the help section of a matlab function
            % that is in the generic format described below. The format
            % consists of a series of sections and sub-sections like this:
            %
            % ------------------------------------------------------------
            % <SECTION NAME i>:
            % <help text for <SECTION NAME i> spanning one or more lines>
            %
            %       -- OR --
            %
            % <SECTION NAME i>:
            % <sub-section1_i>: <help text for sub-section1_i spanning one or more lines>
            % <sub-section2_i>: <help text for sub-section2_i spanning one or more lines>
            %  . . . . . .
            % <sub-sectionN_i>: <help text for sub-sectionN_i spanning one or more lines>
            %
            % ------------------------------------------------------------
            %
            % For <SECTION NAME i> to be considered a section (as opposed to sub-section)
            % it must be in all uppercase, must be followed by a ':', and must not have
            % any other text on that line. A sub-section is any string without white
            % spaces preceding the first occurrence of ':', ' - ' or '--' on a line line.
            %
            % In addition to the generic parsing of sections, this class also provides
            % methods for accessing specific 'designated' sections should they exist in
            % the help comments. A section is defined as designated in this class if it
            % provides public Get* methods to access it. Current designated sections are
            %
            %       Designated Sections    Associated methods
            %       --------------------------------------
            %       'DESCRIPTION:'          % GetDescr()
            %       'INPUT(S):'             % GetParamDescr()
            %       'USAGE OPTIONS:'        % GetUsageOptions()
            %       'PARAMETERS:'           % GetParamUsage()
            %
            % New methods can be added for any new designated sections added to this
            % class in the future
            %
            % Example from Homer3: formal description of user functions
            % help section syntax:
            % --------------------------------------------------------------
            % SYNTAX:
            % [r1,...,rN] = <funcname>(a1,...,aM,p1,...,pL)
            %
            % UI NAME:
            % <User Interface Function Name>
            %
            % DESCRIPTION:
            % <General function description>
            %
            % INPUT:
            % a1: <Description of a1>
            %    . . . . . . . . . .
            % aM: <Description of am>
            % p1: <Description of am>
            %    . . . . . . . . . .
            % pL: <Description of am>
            %
            % OUPUT:
            % r1: <Description of r1>
            %    . . . . . . . . . .
            % rN: <Description of rN>
            %
            % USAGE OPTIONS:
            % <User-friendly name for option 1>: [r11,...,r1N] = <funcname>(a11,...,a1M,p1,...,pL)
            %    . . . . . . . . . .
            % <User-friendly name for option K>: [rK1,...,rKN] = <funcname>(aK1,...,aKM,p1,...,pL)
            %
            % PARAMETERS:
            % p1: [v11,...,v1J]
            %    . . . . . . . . . .
            % pL: [vL1,...,vLJ]
            %
            if isempty(obj.helpstr)
                return;
            end

            sect = obj.FindSections();
            for ii=1:length(sect)
                obj.AddSection(sect{ii});
            end

            % Find the lines in the help that belong to each help section
            obj.FindSectionLines();

            % Remove leading and trailing blank lines from each section
            obj.RemoveBlankLines();

            % Now that we have the lines associated with each help section,
            % assign the text from these lines to corresponding help sections.
            obj.AssignSectionText();

            % We have parsed all sections: now parse sub-sections for each
            % section
            sect = propnames(obj.sections);
            for ii=1:length(sect)
                obj.ParseSubSections(sect{ii});
            end

        end

        % -------------------------------------------------------------
        function sect = FindSections(obj)
            kk=1;
            sect = cell(length(obj.helpstr),1);
            for ii=1:length(obj.helpstr)
                if obj.IsSectionName(obj.helpstr{ii})
                    k = find(obj.helpstr{ii}==':');
                    sect{kk,:} = strtrim(obj.helpstr{ii}(1:k-1));
                    kk=kk+1;
                end
            end
            sect(kk:end) = [];
        end

        % -------------------------------------------------------------
        function FindSectionLines(obj)
            fields = propnames(obj.sections);
            for iLine=1:length(obj.helpstr)
                for jj=1:length(fields)
                    secname = eval( sprintf('obj.sections.%s.name', fields{jj}) );
                    if includes(obj.helpstr{iLine}, secname)
                        EndPrevSection(obj, iLine);
                        eval( sprintf('obj.sections.%s.lines(1) = iLine+1;', fields{jj}) );
                    end
                end
            end
            EndPrevSection(obj, iLine+1);
        end

        % -------------------------------------------------------------
        function AssignSectionText(obj)
            fields = propnames(obj.sections);
            for jj=1:length(fields)
                lines = eval( sprintf('obj.sections.%s.lines(1):obj.sections.%s.lines(2)', fields{jj}, fields{jj}) );
                for iLine = lines
                    if iLine < 1 || isempty(obj.helpstr{iLine})
                        continue;
                    end
                    eval( sprintf('obj.sections.%s.str = sprintf(''%%s%%s\\n'', obj.sections.%s.str, obj.helpstr{iLine});', ...
                        fields{jj}, fields{jj}) );
                end
            end
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Public methods for accessing designated sections
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % -------------------------------------------------------------
        function str = GetUiname(obj)
            str = '';
            if isproperty(obj.sections,'uiname')
                str = strtrim(obj.sections.uiname.str);
            end
        end

        % -------------------------------------------------------------
        function str = GetDescr(obj)
            str = '';
            if isproperty(obj.sections,'description')
                str = obj.sections.description.str;
            end
        end

        % -------------------------------------------------------------
        function str = GetParamDescr(obj, param)
            str = '';
            subsections = [];
            if isproperty(obj.sections,'input')
                subsections = obj.sections.input.subsections;
            end
            if isproperty(obj.sections,'inputs')
                subsections = obj.sections.inputs.subsections;
            end
            if isempty(subsections)
                return;
            end

            % if the input is a number index into the inputs arrays, then we
            % don't need to seach for the param name
            if iswholenum(param)
                if param <= length(subsections)
                    str = subsections(param).str;
                end
                return;
            end

            % Search for the parameter with name <param>
            for ii=1:length(subsections)
                if strcmp(param, subsections(ii).name)
                    str = subsections(ii).str;
                end
            end
        end

        % -------------------------------------------------------------
        function [usage, friendlyname] = GetUsageOptions(obj)
            usage = {};
            friendlyname = {};
            subsections = [];
            if isproperty(obj.sections,'usageoptions')
                subsections = obj.sections.usageoptions.subsections;
            end
            if isempty(subsections)
                return;
            end
            usage        = cell(length(subsections),1);
            friendlyname = cell(length(subsections),1);
            for ii=1:length(subsections)
                friendlyname{ii}  = strtrim(subsections(ii).name);
                usage{ii}         = strtrim(subsections(ii).str);
            end
        end

        % -------------------------------------------------------------
        function [paramname, valformat] = GetParamUsage(obj)
            paramname = {};
            valformat = {};
            subsections = [];
            if isproperty(obj.sections,'parameters')
                subsections = obj.sections.parameters.subsections;
            end
            if isempty(subsections)
                return;
            end
            paramname = cell(length(subsections),1);
            valformat = cell(length(subsections),1);
            for ii=1:length(subsections)
                paramname{ii} = strtrim(subsections(ii).name);
                valformat{ii} = strtrim(subsections(ii).str);
            end
        end

        % -------------------------------------------------------------
        function helpstr = GetStr(obj)
            helpstr = '';
            fields = propnames(obj.sections);
            for ii=1:length(fields)
                name = eval( sprintf('obj.sections.%s.name', fields{ii}) );
                str = eval( sprintf('obj.sections.%s.str', fields{ii}) );
                helpstr = sprintf('%s%s:\n', helpstr, name);
                helpstr = sprintf('%s%s\n', helpstr, str);
            end
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  Methods for parsing sub-sections
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % -------------------------------------------------------------
        function ParseSubSections(obj, section)
            obj.FindSubSectionLines(section);
            obj.AssignSubSectionText(section);
        end

        % -------------------------------------------------------------
        function FindSubSectionLines(obj, section)
            strs = str2cell_fast(eval(sprintf('obj.sections.%s.str', section)), [], 'keepblanks');
            kk=0;
            subsect = [];
            for iLine=1:length(strs)
                name = obj.GetSubSectionName(section, iLine);
                if ~isempty(name)
                    kk=kk+1;

                    % Create new sub-section
                    subsect(kk).name = strtrim(name);
                    subsect(kk).lines = [0,0];
                    subsect(kk).str = '';

                    % Assign starting line to current subsection
                    subsect(kk).lines(1) = iLine;

                    % See if there's sub-section preceding current one. If
                    % so, end it
                    if kk>1 && subsect(kk-1).lines(2)==0
                        subsect(kk-1).lines(2) = iLine-1;
                    end
                end
            end
            if kk>0
                if subsect(kk).lines(2)==0
                    subsect(kk).lines(2) = iLine;
                end
            end
            eval( sprintf('obj.sections.%s.subsections = subsect;', section) );
        end

        % -------------------------------------------------------------
        function [name, k] = GetSubSectionName(obj, section, iLine)
            name = '';
            k = [];
            strs = str2cell_fast(eval(sprintf('obj.sections.%s.str', section)), [], 'keepblanks');

            % Rule 1: Valid section is must end with a ':',' - ', or '--'
            k1 = find(strs{iLine}==':');
            k2 = strfind(strs{iLine}, ' - ');
            k3 = strfind(strs{iLine}, '--');
            d=1;
            if ~isempty(k1)
                k = k1(1);
            elseif ~isempty(k2)
                k = k2(1)+2; d = 2;
            elseif ~isempty(k3)
                k = k3(1)+1;
            end
            if isempty(k)
                return;
            end
            temp = strtrim(strs{iLine}(1:(k-d)));

            %%%% Check to make sure string preceding ':' is a valid section name.

            % Rule 2: If string preceding ':' (ie, potential section name) has spaces then it's
            % not a section name.
            if ~isempty(find(temp==' '))
                return;
            end

            % Rule 3: If string preceding ':' is indented greater than 2 spaces then it's not
            % a section name.
            indent_size = strfind(strs{iLine}(1:(k-d)), temp)-1;
            if isempty(indent_size) || indent_size > 2
                return;
            end

            name = strtrim(strs{iLine}(1:(k-d)));
            k = k(1);
        end

        % -------------------------------------------------------------
        function AssignSubSectionText(obj, section)
            strs = str2cell_fast(eval(sprintf('obj.sections.%s.str', section)), [], 'keepblanks');
            subsect = eval(sprintf('obj.sections.%s.subsections', section));
            for ii=1:length(subsect)
                lines = subsect(ii).lines(1):subsect(ii).lines(2);
                for iLine = lines
                    % The help part corresponding to the subsection name (eg, SD: or trange:)
                    % starts after the colon (:) or double dash (--)
                    [~, k] = obj.GetSubSectionName(section, iLine);
                    if ~isempty(k)
                        strs{iLine} = strtrim(strs{iLine}(k+1:end));
                    end
                    subsect(ii).str = sprintf('%s%s\n', subsect(ii).str , strs{iLine});
                end
            end
            eval( sprintf('obj.sections.%s.subsections = subsect;', section) );
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods (Access = 'private')

        % -----------------------------------------------------------------
        function InitLineIdxs(obj)
            fields = propnames(obj.sections);
            for ii=1:length(fields)
                eval( sprintf('obj.sections.%s.lines = [0 0];', fields{ii}) );
            end
        end

        % -----------------------------------------------------------------
        function RemoveBlankLines(obj)
            % Identify lines in each section where actual text begins and
            % ends and remove the leading and trailing blank lines
            fields = propnames(obj.sections);
            for jj=1:length(fields)
                lines = eval( sprintf('obj.sections.%s.lines(1):obj.sections.%s.lines(2)', fields{jj}, fields{jj}) );

                % If there a zeroes in lines then section is invalid or
                % empty
                if ismember(0, lines)
                    continue;
                end

                % Find first non-blank line of text
                for iLine = lines
                    if ~isblankline(obj.helpstr{iLine})
                        eval( sprintf('obj.sections.%s.lines(1) = iLine;', fields{jj}) );
                        break;
                    end
                end

                % Find last non-blank line of text
                for iLine = fliplr(lines)
                    % Find first non-blank line of text
                    if ~isblankline(obj.helpstr{iLine})
                        eval( sprintf('obj.sections.%s.lines(2) = iLine;', fields{jj}) );
                        break;
                    end
                end
            end
        end

        % -----------------------------------------------------------------
        function EndPrevSection(obj, iLine) %#ok<INUSD>
            fields = propnames(obj.sections);
            for ii=1:length(fields)
                if eval( sprintf('obj.sections.%s.lines(1)>0 && obj.sections.%s.lines(2)==0', fields{ii}, fields{ii}) )
                    eval( sprintf('obj.sections.%s.lines(2) = iLine-1;', fields{ii}) );
                    break;
                end
            end
        end

        % -----------------------------------------------------------------
        function b = IsSectionName(obj, str)
            b = false;
            if isblankline(str)
                return;
            end
            str = strtrim(str);
            if ~isuppercase(str)
                return;
            end
            if str(end)~=':'
                return;
            end
            b = true;
        end

        % ----------------------------------------------------------------------------------
        function Copy(obj, obj2)
            obj.funcname = obj2.funcname;
            obj.helpstr = obj2.helpstr;
            obj.sections = obj2.sections;
        end

        % ----------------------------------------------------------------------------------
        function b = IsEmpty(obj)
            b = false;
            if isempty(obj.funcname)
                b=true;
            end
            if isempty(obj.helpstr)
                b=true;
            end
            if isempty(obj.sections)
                b=true;
            end
        end

        % ---------------------------------------------------------------------------------
        function Help(obj)
            funcname = '';
            if ~isdeployed()
                funcname = [obj.funcname, '.m'];
            else
                for ii = 1:length(obj.userfuncdir)
                    if exist([obj.userfuncdir{ii}, obj.funcname, '.m'], 'file')
                        funcname = [obj.userfuncdir{ii}, obj.funcname, '.m'];
                    end
                end
            end
            if isempty(funcname)
                return
            end
            obj.helpstr = str2cell_fast(help_local(funcname), [], 'keepblanks');
        end

    end

    methods (Static)

        % ----------------------------------------------------------------------------------
        function out = UserFuncDirs(arg)
            persistent userfuncdir;
            if isempty(userfuncdir)
                userfuncdir = FindUserFuncDir(arg);
            end
            out = userfuncdir;
        end

    end

end

author: BUNPC


text:
function variance = english_word_length_variance ( )

%*****************************************************************************80
%
%% ENGLISH_WORD_LENGTH_VARIANCE: variance of the English Word Length PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 July 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Henry Kucera, Winthrop Francis,
%    Computational Analysis of Present-Day American English,
%    Brown University Press, 1967.
%
%  Parameters:
%
%    Output, real VARIANCE, the variance of the PDF.
%
  word_length_max = 27;

  pdf_vec = [ ...
    0.03160, ...
    0.16975, ...
    0.21192, ...
    0.15678, ...
    0.10852, ...
    0.08524, ...
    0.07724, ...
    0.05623, ...
    0.04032, ...
    0.02766, ...
    0.01582, ...
    0.00917, ...
    0.00483, ...
    0.00262, ...
    0.00099, ...
    0.00050, ...
    0.00027, ...
    0.00022, ...
    0.00011, ...
    0.00006, ...
    0.00005, ...
    0.00002, ...
    0.00001, ...
    0.00001, ...
    0.00001, ...
    0.00001, ...
    0.00001 ];
  pdf_sum = 0.99997;

  mean = 0.0;
  for j = 1 : word_length_max
    mean = mean + j * pdf_vec(j);
  end

  mean = mean / pdf_sum;

  variance = 0.0;
  for j = 1 : word_length_max
    variance = variance + pdf_vec(j) * ( j - mean )^2;
  end

  variance = variance / pdf_sum;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [F,feasible] = filter_norm_2(all_f,all_c_w,all_c_x,all_Q_xw,x,Zmodel,allw,X,Q_xx,VariableType)

feasible = 1;
% As a first step, we figure out the radius
r = Zmodel.r;
center = Zmodel.center;
F = ([]);
ci_basis = all_c_w';
lastBici = [];
lastusedrows = [];
% (bi' + (Bi*w)')*x + (ci'*w + di).
for i = 1:length(all_f)
    Bi = 2*all_Q_xw(:,length(x)*(i-1)+1:length(x)*i)';
    bi = all_c_x(length(x)*(i-1)+1:length(x)*i);
    if (nnz(ci_basis(:,i))==0) & nnz(Bi)==0
        F = F + (X(i)>=0);
        % This constraint row is constant
    else
        ci = ci_basis(:,i);
        di = all_f(i);
        used = find(full(any([full(Bi') ci],2)));
        %        ci = ci(used);
        %        Bi = Bi(:,used);
        % Shift |w-center|, wtilde = w-center i.e. w=wtilde+center
        di = di + ci'*center;
        bi = bi + Bi*center;
        if isequal(abs(lastBici),abs([Bi' ci]))
             F = F + (x'*Q_xx{i}*x+bi'*x + di - r*s >= 0);
      %       F = F + (x'*Q_xx{i}*x+bi'*x + di - r*norm(full(Bi(used_rows,used)))*s >= 0);
        else
            s = sdpvar(1,1);
            if length(used)==1
                if nnz(Bi(:,used))==0
                    s = norm(full(ci(used)));
                    F = F + (x'*Q_xx{i}*x+bi'*x + di - r*s >= 0);
                else
                    F = F + (x'*Q_xx{i}*x+bi'*x + di - r*s >= 0) + (-s<=Bi(:,used)'*x+ci(used)<=s);
                end
                used_rows = 1:size(Bi,1);
                lastBici = [Bi' ci];
            else
                used_rows = find(any(full(Bi(:,used)')));
                if length(used_rows)==1 & nnz(ci(used))==0
                    % Special case norm(double vector*scalar)
                    if isequal(used_rows,lastusedrows)
                        F = F + (x'*Q_xx{i}*x+bi'*x + di - r*norm(full(Bi(used_rows,used)))*lasts >= 0)  ;
                    else
                        F = F + (x'*Q_xx{i}*x+bi'*x + di - r*norm(full(Bi(used_rows,used)))*s >= 0) + (-s <= x(used_rows)<=s);
                        lastBici = [Bi' ci];
                        lastusedrows = used_rows;
                        lasts = s;
                    end
                else
                    if nnz(Bi)==0
                        s = norm(full(ci(used)));
                        F = F + (x'*Q_xx{i}*x+bi'*x + di - r*s >= 0) ;
                    else
                        F = F + (x'*Q_xx{i}*x+bi'*x + di - r*s >= 0) + (cone(Bi(:,used)'*x+ci(used),s));
                    end
                    lastBici = [Bi' ci];
                end
            end
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function theta = circle3dPosition(point, circle)
%CIRCLE3DPOSITION Return the angular position of a point on a 3D circle
%
%   POS = circle3dPosition(POINT, CIRCLE)
%   Returns angular position of point on the circle, in degrees, between 0
%   and 360.
%   with POINT: [xp yp zp]
%   and CIRCLE: [X0 Y0 Z0 R THETA PHI] or [X0 Y0 Z0 R THETA PHI PSI]
%   (THETA being the colatitude, and PHI the azimut)
%
%   See also:
%   circles3d, circle3dOrigin, circle3dPoint
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/02/2005
%

%   HISTORY
%   27/06/2007: change 3D angle convention
%   2011-06-21 use degrees for angles

% get center and radius
xc = circle(:,1);
yc = circle(:,2);
zc = circle(:,3);

% get angle of normal
theta   = circle(:,5);
phi     = circle(:,6);

% find origin of the circle
ori     = circle3dOrigin(circle);

% normal vector of the supporting plane (cartesian coords)
vn      = sph2cart2d([theta phi]);

% create plane containing the circle
plane   = createPlane([xc yc zc], vn);

% find position of point on the circle plane
pp0     = planePosition(ori,    plane);
pp      = planePosition(point,  plane);

% compute angles in the planes
theta0  = mod(atan2(pp0(:,2), pp0(:,1)) + 2*pi, 2*pi);
theta   = mod(atan2(pp(:,2), pp(:,1)) + 2*pi - theta0, 2*pi);

% convert to degrees
theta = theta * 180 / pi;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [Gs, op, nodes] = mk_nbrs_of_digraph(G0)
% MK_NBRS_OF_DIGRAPH Make all digraphs that differ from G0 by a single edge deletion, addition or reversal
% [Gs, op, nodes] = mk_nbrs_of_digraph(G0)
%
% Gs(:,:,i) is the i'th neighbor
% op{i} = 'add', 'del', or 'rev' is the operation used to create the i'th neighbor.
% nodes(i,1:2) are the head and tail of the operated-on arc.

debug = 0; % the vectorized version is about 3 to 10 times faster

n = length(G0);
[I,J] = find(G0); % I(k), J(k) is the k'th edge
E = length(I); % num edges present in G0

% SINGLE EDGE DELETIONS

Grep = repmat(G0(:), 1, E); % each column is a copy of G0
% edge_ndx(k) is the scalar location of the k'th edge
edge_ndx = find(G0);
% edge_ndx = subv2ind([n n], [I J]); % equivalent
% We set (ndx(k), k) to 0 for k=1:E in Grep
ndx = subv2ind(size(Grep), [edge_ndx(:) (1:E)']);
G1 = Grep;
G1(ndx) = 0;
Gdel = reshape(G1, [n n E]);

% if debug
% % Non-vectorized version
% ctr = 1;
% for e=1:E
%   i = I(e); j = J(e);
%   Gdel2(:,:,ctr) = G0;
%   Gdel2(i,j,ctr) = 0;
%   ctr = ctr + 1;
% end
% assert(isequal(Gdel, Gdel2));
% end

% SINGLE EDGE REVERSALS

% rev_edge_ndx(k) is the scalar location of the k'th reversed edge
%rev_edge_ndx = find(G0'); % different order to edge_ndx, which is bad
rev_edge_ndx = subv2ind([n n], [J I]);
% We set (rev_edge_ndx(k), k) to 1 for k=1:E in G1
% We have already deleted i->j in the previous step
ndx = subv2ind(size(Grep), [rev_edge_ndx(:) (1:E)']);
G1(ndx) = 1;
Grev = reshape(G1, [n n E]);

% if debug
% % Non-vectorized version
% ctr = 1;
% for e=1:E
%   i = I(e); j = J(e);
%   Grev2(:,:,ctr) = G0;
%   Grev2(i,j,ctr) = 0;
%   Grev2(j,i,ctr) = 1;
%   ctr = ctr + 1;
% end
% assert(isequal(Grev, Grev2));
% end

% SINGLE EDGE ADDITIONS

Gbar = ~G0; % Gbar(i,j)=1 iff there is no i->j edge in G0
Gbar = setdiag(Gbar, 0); % turn off self loops
[Ibar,Jbar] = find(Gbar);

bar_edge_ndx = find(Gbar);
Ebar = length(Ibar); % num edges present in Gbar
Grep = repmat(G0(:), 1, Ebar); % each column is a copy of G0
ndx = subv2ind(size(Grep), [bar_edge_ndx(:) (1:Ebar)']);
Grep(ndx) = 1;
Gadd = reshape(Grep, [n n Ebar]);

% if debug
% % Non-vectorized version
% ctr = 1;
% for e=1:length(Ibar)
%   i = Ibar(e); j = Jbar(e);
%   Gadd2(:,:,ctr) = G0;
%   Gadd2(i,j,ctr) = 1;
%   ctr = ctr + 1;
% end
% assert(isequal(Gadd, Gadd2));
% end

Gs = cat(3, Gdel, Grev, Gadd);

nodes = [I J;
	 I J;
	 Ibar Jbar];

op = cell(1, E+E+Ebar);
op(1:E) = {'del'};
op(E+1:2*E) = {'rev'};
op(2*E+1:end) = {'add'};

% numeric output:
% op(i) = 1, 2, or 3, if the i'th neighbor was created by adding, deleting or reversing an arc.

ADD = 1;
DEL = 2;
REV = 3;

%op = [repmat(DEL, 1, E) repmat(REV, 1, E) repmat(ADD, 1, Ebar)];

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function table2 = i4mat_border_add ( m, n, table )

%*****************************************************************************80
%
%% I4MAT_BORDER_ADD adds a "border" to an I4MAT.
%
%  Discussion:
%
%    We suppose the input data gives values of a quantity on nodes
%    in the interior of a 2D grid, and we wish to create a new table
%    with additional positions for the nodes that would be on the
%    border of the 2D grid.
%
%                  0 0 0 0 0 0
%      * * * *     0 * * * * 0
%      * * * * --> 0 * * * * 0
%      * * * *     0 * * * * 0
%                  0 0 0 0 0 0
%
%    The illustration suggests the situation in which a 3 by 4 array
%    is input, and a 5 by 6 array is to be output.
%
%    The old data is shifted to its correct positions in the new array.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the table data.
%
%    Output, integer TABLE2(M+2,N+2), the augmented table data.
%
  table2 = zeros(m+2,n+2);

  table2(1,1:n+2) = 0;
  table2(m+2,1:n+2) = 0;
  table2(2:m+1,1) = 0;
  table2(2:m+1,n+2) = 0;

  table2(2:m+1,2:n+1) = round ( table(1:m,1:n) );

  return
end


text:
function [model, B, elapse] = LSH_learn(A, maxbits)
%   This is a function of LSH (Locality Sensitive Hashing) learning.
%
%	Usage:
%	[model, B,elapse] = LSH_learn(A, maxbits)
%
%	      A: Rows of vectors of data points. Each row is sample point
%   maxbits: Code length
%
%     model: Used for encoding a test sample point.
%	      B: The binary code of the input data A. Each row is sample point
%    elapse: The coding time (training time).
%
%
%
%   Reference:
%
%   Moses S. Charikar: Similarity estimation techniques from rounding
%   algorithms. Proceedings of the thiry-fourth annual ACM symposium on
%   Theory of computing, 2002.
%
%
%   version 2.0 --Nov/2016
%   version 1.0 --Jan/2010
%
%   Written by  Yue Lin (linyue29@gmail.com)
%               Deng Cai (dengcai AT gmail DOT com)
%

tmp_T = tic;

[~,Nfeatures] = size(A);
k = maxbits;

U = normrnd(0, 1, Nfeatures, k);
Z = A * U;

B = (Z > 0);
model.U = U;

elapse = toc(tmp_T);
end

author: ZJULearning

repo: MatlabFunc


text:
classdef prtKernelSet < prtKernel
    % prtKernelSet  A set of prtKernel Object
    %
    %   KC = prtKernelSet(KERNOBJ) creates a prtKernelSet object consisting
    %   of one prtKenrel KERNOBJ.
    %
    %   KC = KERNOBJ1 & KERNOBJ2 creates a prtKernelSet object consisting of
    %   the prtKernel objects KERNOBJ1, and KERNOBJ2. Additional prtKernel
    %   object can be included by appending more objects using the AND(&)
    %   operator.
    %
    %   When the TRAIN method is called on a prtKernetSet object, the
    %   train method of all member prtKernel objects is invoked.
    %
    %   When the RUN method is called on a prtKernelSet object is called,
    %   the RUN method of all member prtKernel objects is invoked, and the
    %   output is the concatenation of the results.
    %
    %   % Example:
    %
    %   ds = prtDataGenBimodal;       % Generate a dataset
    %   k1 = prtKernelRbf;            % Create a prtKernel object with
    %                                 % default value of sigma
    %   k2 = prtKernelRbf('sigma',2); % Create a prtKernel object with
    %                                 % specified value of sigma
    %
    %   kSet = k1&k2;                 % Create a prtKernel object using the
    %                                 % AND operator
    %   kSet =  kSet.train(ds);       % Train
    %
    %   g = kSet.run(ds);             % Evaluate
    %
    %   imagesc(g.getObservations);   %Plot the results
    %
    %   See also: prtKernel,prtKernelSet, prtKernelDc, prtKernelDirect,
    %   prtKernelHyperbolicTangent, prtKernelPolynomial,
    %   prtKernelRbfNdimensionScale

    properties (SetAccess = private)
        name = 'Kernel Set'; % Kernel Set
        nameAbbreviation = 'KernelSet'; % KernelSet
    end

    properties
        kernelCell % The cell array of prtKenrel objects
    end

    methods (Access = protected, Hidden = true)

        function Obj = trainAction(Obj,ds)
            for i = 1:length(Obj.kernelCell)
                Obj.kernelCell{i} = Obj.kernelCell{i}.train(ds);
            end
            Obj.isTrained = true;
        end

        function dsOut = runAction(Obj,ds)

            for i = 1:length(Obj.kernelCell)
                if i == 1
                    dsOut = Obj.kernelCell{i}.run(ds);
                else
                    %Note: use getObservations here to make sure that
                    %there's no confusion when cat'ing features from
                    %prtDataSetClass and prtDataSetRegress'es
                    dsOut = dsOut.catFeatures(getObservations(Obj.kernelCell{i}.run(ds)));
                end
            end
        end
    end
    methods
        function Obj = prtKernelSet(varargin)
            %Obj = prtKernelSet(varargin)
            %   prtKernelSet(kernel1,kernel2,...)

            if nargin == 0
                error('prt:prtKernelSet:tooFewInputs','prtKernelSet requires at least one prtKernel object as input');
            end
            c = 1;
            tempKernelCell = {};
            for i = 1:length(varargin)
                if isa(varargin{i},'prtKernelSet')
                    tempCell = varargin{i}.getKernelCell;
                    tempKernelCell(c:c+length(tempCell)-1) = tempCell;
                    c = c + length(tempCell);
                elseif isa(varargin{i},'prtKernel')
                    tempKernelCell{c} = varargin{i}; %#ok<AGROW>  Impossible to tell how big this will get...
                    c = c + 1;
                end
            end
            Obj.kernelCell = tempKernelCell;
        end

        function Obj = set.kernelCell(Obj,aCell)

            %Check is cell:
            if ~isa(aCell,'cell')
                error('prtKernel:kernelCell','prtKernel''s kernelCell must be a cell array');
            end
            %Check right size:
            if ~isvector(aCell)
                error('prtKernel:kernelCell','prtKernel''s kernelCell must be a vector cell array');
            end
            %             if length(aCell) ~= size(Obj.connectivityMatrix,1)-2
            %                 error('prtKernel:kernelCell','Attempt to change a prtKernel''s kernelCell''s size.  kernelCell must be a vector cell array of length(size(Obj.connectivityMatrix,1)-2)');
            %             end

            %Check all are prtKernels:
            if ~all(cellfun(@(c)isa(c,'prtKernel'),aCell))
                error('prtKernel:kernelCell','kernelCell must be a vector cell array of prtKernels')
            end

            %Set the internal action cell correctly
            Obj.kernelCell = aCell;
        end

    end

    methods (Hidden = true)
        function [nDimensions,nDimensionsArray] = nDimensions(Obj)

            nDimensionsArray = nan(length(Obj.kernelCell),1);
            for i = 1:length(Obj.kernelCell)
                nDimensionsArray(i) = Obj.kernelCell{i}.nDimensions;
            end
            nDimensions = sum(nDimensionsArray);
        end

        function Obj = retainKernelDimensions(Obj,keepLogical)

            if ~islogical(keepLogical)
                temp = false(1,Obj.nDimensions);
                temp(keepLogical) = true;
                keepLogical = temp;
            end

            start = 1;
            for i = 1:length(Obj.kernelCell)
                nDimensions = Obj.kernelCell{i}.nDimensions;
                Obj.kernelCell{i} = Obj.kernelCell{i}.retainKernelDimensions(keepLogical(start:start+nDimensions-1));
                start = start+nDimensions;
            end
        end

    end

    methods (Hidden = true)
        function h = plot(Obj)
            %Plot each kernel:
            h = zeros(length(Obj.kernelCell),1);
            for i = 1:length(Obj.kernelCell)
                hT = Obj.kernelCell{i}.plot;
                if ~isempty(hT)
                    h(i) = hT;
                end
            end
        end
    end

    methods (Hidden = true)
        function kernelCell = getKernelCell(Obj)
            kernelCell = Obj.kernelCell;
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% DECIMATE_LIBIGL Decimate a closed manifold mesh (V,F)
%
% [W,G] = decimate_libigl(V,F,ratio)
% [W,G,J,I] = decimate_libigl(V,F,ratio,'ParameterName',ParameterValue, ...)
%
% Inputs:
%   V  #V by 3 list of vertex positions
%   F  #F by 3 list of triangle indices into V
%   ratio   either a 1<number<#F  of max faces, or a 0<ratio<1 to be multiplied
%     against #F to get max faces in output
%   Optional:
%     'Method' followed by one of:
%        {'naive'}  simply collapse small edges and place vertices at midpoint
%        'qslim'  Quadric error metric
%        'progressive-hulls'  Progressive hulls
% Outputs:
%   W  #W by 3 list of vertex positions
%   G  #G by 3 list of triangle indices into W
%   J  #G list of indices into F of birth face
%   I  #U list of indices into V of birth vertices
%

author: alecjacobson

repo: gptoolbox


text:
function B = mirror_shrink(A)
    [m,n] = size(A);
    yi = 2:m-1;
    xi = 2:n-1;
    B = A(yi,xi);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
classdef LiftForceModel < AbstractForceModel
    %LiftForceModel Summary of this class goes here
    %   Detailed explanation goes here

    properties

    end

    methods
        function obj = LiftForceModel()

        end

        function [forceVect, tankMdots, ecStgDots] = getForce(obj, ut, rVect, vVect, mass, bodyInfo, aero, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, attState)
            if(norm(rVect) - (bodyInfo.radius + bodyInfo.atmohgt) > 0)
                forceVect = [0;0;0];
            else
                forceVect = getLiftForce(bodyInfo, ut, rVect, vVect, aero, mass, attState);
            end

            tankMdots = [];
            ecStgDots = [];
        end
    end
end

function forceVect = getLiftForce(bodyInfo, ut, rVectECI, vVectECI, aero, mass, attState)
    arguments
        bodyInfo(1,1) KSPTOT_BodyInfo
        ut(1,1) double
        rVectECI(3,1) double
        vVectECI(3,1) double
        aero(1,1) LaunchVehicleAeroState
        mass(1,1) double
        attState(1,1) LaunchVehicleAttitudeState
    end

    rVectECI = reshape(rVectECI,3,1);
    vVectECI = reshape(vVectECI,3,1);

    altitude = norm(rVectECI) - bodyInfo.radius;

    if(altitude <= bodyInfo.atmohgt && altitude >= 0)
        [lat, long, ~, ~, ~, ~, ~, vVectECEF] = getLatLongAltFromInertialVect(ut, rVectECI, bodyInfo, vVectECI);
        [density, pressureKPA, ~] = getAtmoDensityAtAltitude(bodyInfo, altitude, lat, ut, long);
    elseif(altitude <= 0)
        density = 0;
    else
        density = 0;
    end

    if(density > 0)
        [ClS, liftUnitVectInertial] = aero.getLiftCoeffAndDir(ut, rVectECI, vVectECI, bodyInfo, mass, altitude, pressureKPA, density, vVectECEF, attState);

        vVectECEFMag = norm(vVectECEF);

        %all forces are returned in units of mT*km/s^2 = (1000000)*N
        FL = (1/2)*density*(vVectECEFMag^2)*ClS; %kg/m^3 * (km^2/s^2) * m^2 = kg/m * km^2/s^2 = kg*(km/m)*km/s^2 = kg*(1000)*km/s^2

        forceVect = FL * liftUnitVectInertial;

%         Cl_level = aero.Cl_0;
%         A = aero.areaLift;
%         bodyLiftVect = normVector(aero.bodyLiftVect);
%
%         vVectEcefMag = norm(vVectECEF);
%
%         body2InertDcm = steeringModel.getBody2InertialDcmAtTime(ut, rVectECI, vVectECI, bodyInfo);
%         [~,angOfAttack,~] = computeAeroAnglesFromBodyAxes(ut, rVectECI, vVectECI, bodyInfo, body2InertDcm(:,1), body2InertDcm(:,2), body2InertDcm(:,3));
%
%         Cl = Cl_level + 2*pi*angOfAttack;
%
%         ClA = Cl * A;
%
%         L = (1/2) * density * (vVectEcefMag^2) * ClA; %kg/m^3 * (km^2/s^2) * m^2 = kg/m * km^2/s^2 = kg*(km/m)*km/s^2 = kg*(1000)*km/s^2
%
%         bodyLiftForce = L * bodyLiftVect;
%
%         liftForce = body2InertDcm * bodyLiftForce;
    else
        forceVect = [0;0;0];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function h = overlayMeanTSeries(view,scan1,scan2);
%  h = overlayMeanTSeries(view,[scan1,scan2]);
%
% Provides an interface for viewing mean tSeries on top of one
% another.Created for the purpose of seeing if movement occurred.
%
% scan1 and scan2 are the two tSeries to compare. Will compute the
% mean (all slices), and overlay them. They default to the first
% and last scan in the current data type.
%
% If either scan1 or scan2 is entered as 0, this will cause
% a dialog to be brought up to get the scan numbers.
%
% written 06/17/04 ras.
if ieNotDefined('scan1')
    scan1 = 1;
end

if ieNotDefined('scan2')
    dt = viewGet(view,'curdt');
    global dataTYPES;
    scan2 = length(dataTYPES(dt).scanParams);
end

slices = 1:numSlices(view);

if scan1==0 | scan2==0
    % get the scan nums with a dialog
    prompt = {'First Scan:','Second Scan:'};
    dlgTitle = 'Overlay Mean functionals from which scans?';
    defaults = {'1' '2'};
    answer = inputdlg(prompt,dlgTitle,1,defaults);
    scan1 = str2num(answer{1});
    scan2 = str2num(answer{2});
end

hbox = msgbox('Getting mean tSeries to overlay...');

mapPath = fullfile(dataDir(view),'meanMap.mat');

if ~exist(mapPath,'file')
    % compute the mean map
    mrGlobals;
    loadSession;
    hI = initHiddenInplane;
    hI = computeMeanMap(hI,0);
    clear hI;
end

load(mapPath,'map');

firstScan = map{scan1};
lastScan = map{scan2};

name1 = sprintf('tSeries Scan %i',scan1);
name2 = sprintf('tSeries Scan %i',scan2);

overlayMaps(view,firstScan,lastScan,1,name1,name2);

close(hbox);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [prevhold] = takehold(ax, mode)
%TAKEHOLD  hold axes and return previous hold status.
%   TAKEHOLD(ax) is equivalent to hold(ax, 'on'), where ax is the handle to
%   an axes object.
%
%   TAKEHOLD(ax, 'local') clears the axes ax if their previous hold state
%   is off and then holds the axes, i.e., hold(ax, 'on').
%   This is useful when a plotting function draws multiple
%   graphics objects, so it needs to hold the plot intermediately, but
%   maintains pre-existing graphics in the plot only if hold(ax 'on') has
%   been issued by another source. This is the same behavior that the plot
%   function exhibits.
%
%   TAKEHOLD(ax, 'add') or TAKEHOLD(ax) does not clear the axes ax, but
%   instead holds it. It adds to the plot, without erasing previous
%   graphics.
%
%   prevhold = TAKEHOLD(ax, ...) returns the previous hold status of the
%   axes object ax. This can then be restored by calling
%
% usage
%   prevhold = TAKEHOLD
%   prevhold = TAKEHOLD(ax)
%   prevhold = TAKEHOLD(ax, 'add')
%   prevhold = TAKEHOLD(ax, 'local')
%
% input
%   ax = axes object handle
%   mode = 'add' | 'local'
%
% output
%   held = 1 if ax already held
%          0 if ax not already held
%
% See also RESTOREHOLD, HOLD, ISHOLD.
%
% File:      takehold.m
% Author:    Ioannis Filippidis, jfilippidis@gmail.com
% Date:      2012.01.22 - 2012.05.10
% Language:  MATLAB R2012a
% Purpose:   hold(ax, 'on') and return current hold status
% Copyright: Ioannis Filippidis, 2012-

% which axes' handle ?
if nargin < 1
    ax = gca;
end

% which mode ?
if nargin < 2
    mode = 'add';
end

prevhold = ishold(ax);

if prevhold == 1
    return
end

if strcmp(mode, 'local')
    disp('heh')
    cla(ax)
end

hold(ax, 'on')

author: Sable

repo: mcbench-benchmarks


text:
function result  = omp(Dict, X, K, epsilon, options)
    if(isobject(Dict))
        Dict = double(Dict);
    end
    if ~isa(Dict, 'double')
        error('Dict must be a double matrix');
    end
    if nargin < 4
        epsilon = 1e-3;
    end
    if nargin < 5
        options = struct;
    end
    if ~isfield(options, 'sparse_output')
        options.sparse_output = 1;
    end
    % Options for least squares
    ls_ls = 0;
    ls_chol = 1;
    if ~isfield(options, 'ls_method')
        options.ls_method = 'chol';
    end
    ls_method = ls_chol;
    if strcmp(options.ls_method,'ls')
        ls_method = ls_ls;
    elseif strcmp(options.ls_method, 'chol')
        ls_method = ls_chol;
    end
    if ~isfield(options, 'verbose')
        options.verbose = 0;
    end
    result = mex_omp_chol(Dict, X, K, epsilon, ...
        options.sparse_output,...
        ls_method, ...
        options.verbose);
end

author: indigits

repo: sparse-plex


text:
function [remoteDoses, planC] = listRemoteDoseArrays(planC, saveFlag)
%"listRemoteDoseArrays"
%   Returns a structure array of dose array variables that are stored
%   remotely, for the purpose of wrapping the remote files in the archive
%   when it is saved, as well as updating the remote files.
%
%   If saveFlag is set to 1 the contents of the data field for each
%   remotely stored variable is written to that variable's file.
%   The variable is then marked as unloaded and returned as part of
%   the modified planC.  By default this is off.
%
%   If no planC is passed, listRemoteDoseArrays operates on the global
%   planC.
%
%JRA 04/26/05
%
%Usage:
%   [remoteDoses, planC] = listRemoteDoseArrays(planC, saveFlag)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('saveFlag')
    saveFlag = 0;
end

if ~exist('planC')
    global planC
end
indexS = planC{end};

remoteDoses = [];
nDoses = length(planC{indexS.dose});

%Iterate over doses, checking for remote dose arrays.
for doseNum=1:nDoses

    %Get the dose array for future checking.
    structure = planC{indexS.dose}(doseNum).doseArray;

    if ~isLocal(structure)

        %Found a remote array, add it to our list.
        remoteDoses = [remoteDoses structure];

        %If the saveFlag is specified, write the data field to the file.
        if saveFlag
            if structure.isLoaded
                structure = saveRemoteVariable(structure);
            end
        end
    end
end
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%% Copyright (C) 2017, 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @deftypefun {@var{y} =} double_to_sym_exact (@var{x})
%% Convert a double value to the equivalent rational sym
%%
%% Private helper function.
%%
%% @end deftypefun

function y = double_to_sym_exact (x)
  if (isnan (x))
    y = pycall_sympy__ ('return S.NaN');
  elseif (isinf (x) && x < 0)
    y = pycall_sympy__ ('return -S.Infinity');
  elseif (isinf (x))
    y = pycall_sympy__ ('return S.Infinity');
  else
    %% Rational will exactly convert from a float
    y = pycall_sympy__ ('return Rational(_ins[0])', x);
  end
end

repo: octsympy


text:
classdef MimEditGrey < MimGuiPlugin
    % MimEditGrey. Gui Plugin for setting paint colour
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimEditGrey is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Grey'
        SelectedText = 'Grey'
        ToolTip = 'Edit with grey label'
        Category = 'Segmentation label'
        Visibility = 'Dataset'
        Mode = 'Edit'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 5
        ButtonHeight = 1
        Icon = 'paint.png'
        IconColour = GemMarkerPoint.DefaultColours{7}
        Location = 37
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.PaintBrushColour = 7;
        end

        function enabled = IsEnabled(gui_app)
             enabled = gui_app.IsDatasetLoaded && gui_app.ImagePanel.OverlayImage.ImageExists && ...
                isequal(gui_app.ImagePanel.SelectedControl, 'Paint');
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.PaintBrushColour == 7;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
classdef PluginConstraint < AbstractConstraint
    %PluginConstraint Summary of this class goes here
    %   Detailed explanation goes here

    properties
        normFact = 1;
        plugin LvdPlugin
        event LaunchVehicleEvent
        eventNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        evalType(1,1) ConstraintEvalTypeEnum = ConstraintEvalTypeEnum.FixedBounds;
        stateCompType(1,1) ConstraintStateComparisonTypeEnum = ConstraintStateComparisonTypeEnum.Equals;
        stateCompEvent LaunchVehicleEvent
        stateCompNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;
    end

    methods
        function obj = PluginConstraint(plugin, event, lb, ub)
            obj.plugin = plugin;
            obj.event = event;
            obj.lb = lb;
            obj.ub = ub;

             obj.id = rand();
        end

        function [lb, ub] = getBounds(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function [c, ceq, value, lwrBnd, uprBnd, type, eventNum, valueStateComp] = evalConstraint(obj, stateLog, celBodyData)
            type = obj.getConstraintType();

            switch obj.eventNode
                case ConstraintStateComparisonNodeEnum.FinalState
                    stateLogEntry = stateLog.getLastStateLogForEvent(obj.event);

                case ConstraintStateComparisonNodeEnum.InitialState
                    stateLogEntry = stateLog.getFirstStateLogForEvent(obj.event);

                otherwise
                    error('Unknown event node.');
            end

            if(not(isempty(obj.frame)))
                frame = obj.frame;
            else
                frame = stateLogEntry.centralBody.getBodyCenteredInertialFrame();
            end

            lvdData = stateLogEntry.lvdData;
            pluginSet = lvdData.plugins;
            value = obj.plugin.executePlugin(stateLogEntry.lvdData, stateLog, obj.event, LvdPluginExecLocEnum.Constraint, [],[],[], pluginSet.userData, stateLogEntry, frame);

            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                switch obj.stateCompNode
                    case ConstraintStateComparisonNodeEnum.FinalState
                        stateLogEntryStateComp = stateLog.getLastStateLogForEvent(obj.stateCompEvent);

                    case ConstraintStateComparisonNodeEnum.InitialState
                        stateLogEntryStateComp = stateLog.getFirstStateLogForEvent(obj.stateCompEvent);

                    otherwise
                        error('Unknown event node.');
                end

                valueStateComp = obj.plugin.executePlugin(lvdData, stateLog, obj.event, LvdPluginExecLocEnum.Constraint, [],[],[], pluginSet.userData, stateLogEntryStateComp, frame);
            else
                valueStateComp = NaN;
            end

            [c, ceq] = obj.computeCAndCeqValues(value, valueStateComp);

            lwrBnd = obj.lb;
            uprBnd = obj.ub;

            eventNum = obj.event.getEventNum();
        end

        function sF = getScaleFactor(obj)
            sF = obj.normFact;
        end

        function setScaleFactor(obj, sF)
            obj.normFact = sF;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesEvent(obj, event)
            tf = obj.event == event;
            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                tf = tf || obj.stateCompEvent == event;
            end
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function tf = usesExtremum(obj, extremum)
            tf = false;
        end

        function tf = usesGroundObj(obj, grdObj)
            tf = false;
        end

        function tf = usesPlugin(obj, plugin)
            tf = obj.plugin == plugin;
        end

        function tf = canUseSparseOutput(obj)
            tf = true;
        end

        function event = getConstraintEvent(obj)
            event = obj.event;
        end

        function type = getConstraintType(obj)
            type = 'Plugin Value';
        end

        function [unit, lbLim, ubLim, usesLbUb, usesCelBody, usesRefSc] = getConstraintStaticDetails(obj)
            unit = ''; %TODO: Need to come up with a way to populate this I think?
            lbLim = -Inf;
            ubLim = Inf;
            usesLbUb = true;
            usesCelBody = false;
            usesRefSc = false;
        end

        function addConstraintTf = openEditConstraintUI(obj, lvdData)
            if(lvdData.plugins.getNumPlugins() >= 1)
%                 addConstraintTf = lvd_EditPluginConstraintGUI(obj, lvdData);

                output = AppDesignerGUIOutput({false});
                lvd_EditPluginConstraintGUI_App(obj, lvdData, output);
                addConstraintTf = output.output{1};
            else
                errordlg('There are currently no plugins in this scenario.  Create at least one new plugin first.');

                addConstraintTf = false;
            end
        end

        function plugin = selectConstraintObj(obj, lvdData)
            [listBoxStr, plugins] = lvdData.plugins.getListboxStr();

            plugin = [];
            if(isempty(plugins))
                warndlg('Cannot create plugin value object: no plugins have been created.  Create a plugin first.','Plugin Value Constraint','modal');
            else
                [Selection,ok] = listdlg('PromptString',{'Select a plugin:'},...
                                'SelectionMode','single',...
                                'Name','Plugins',...
                                'ListString',listBoxStr);

                if(ok == 0)
                    plugin = [];
                else
                    plugin = plugins(Selection);
                end
            end
        end

        function useObjFcn = setupForUseAsObjectiveFcn(obj,lvdData)
            pluginSel = obj.selectConstraintObj(lvdData);

            if(not(isempty(pluginSel)))
                obj.plugin = pluginSel;
                useObjFcn = true;
            else
                useObjFcn = false;
            end
        end
    end

    methods(Static)
        function constraint = getDefaultConstraint(~, ~)
            constraint = PluginConstraint(LvdPlugin.empty(1,0), LaunchVehicleEvent.empty(1,0), 0, 0);
        end
    end
end


text:
% PLUGUIN_UIUPDATE - support function for plugin_menu

% Copyright (C) 2019 Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function plugin_uiupdate(fig)

tmpobj = get(fig, 'userdata');
listboxVal = get(findobj(fig, 'tag', 'pluginlist'), 'value');

allPlugins = tmpobj.allplugins;
selectedPlugins = tmpobj.selectedplugins;
selection       = tmpobj.selection;
if length(listboxVal) == 1
    lastSelect = listboxVal;
else
    lastSelect = setdiff(listboxVal, selection);

    % check if selection is consistent
    if ~all([selectedPlugins([ listboxVal ]).installed]) && ~all(~[selectedPlugins([ listboxVal ]).installed] )
        listboxVal = lastSelect(1);
        set(findobj(fig, 'tag', 'pluginlist'), 'value', listboxVal);
    end
end
if isempty(lastSelect), lastSelect = listboxVal(1); end
if length(lastSelect) > 1, lastSelect = lastSelect(1); end

tmpstr = selectedPlugins(lastSelect).description;

% set text
hDescription = findobj(fig, 'tag', 'description');
hTags        = findobj(fig, 'tag', 'tags');
hStatus      = findobj(fig, 'tag', 'status');
hSize        = findobj(fig, 'tag', 'size');
tmpstr = textwrap(hDescription, {tmpstr});
set(hDescription, 'string', tmpstr);
set(hTags       , 'string', selectedPlugins(lastSelect).rawtags );
set(hSize       , 'string', [num2str(ceil(selectedPlugins(lastSelect).size/100)/10) ' MB'] );

% installed status
if selectedPlugins(lastSelect).installed
    if selectedPlugins(lastSelect).installorupdate
        set(hStatus, 'string', ['version ' selectedPlugins(lastSelect).currentversion ' installed - an update is available'] );
    else
        set(hStatus , 'string', 'installed' );
    end
else
    set(hStatus, 'string', 'not installed');
end

% set buttons
set(findobj(fig, 'tag', 'install'), 'enable', 'on');
if selectedPlugins(lastSelect).installed
    set(findobj(fig, 'tag', 'rmbut'), 'enable', 'on');
else
    set(findobj(fig, 'tag', 'rmbut'), 'enable', 'off');
end
if ~isempty(selectedPlugins(lastSelect).webdoc)
    set(findobj(fig, 'tag', 'documentation'), 'enable', 'on', 'callback', ['web(''' selectedPlugins(lastSelect).webdoc ''', ''-browser'');' ]);
else
    set(findobj(fig, 'tag', 'documentation'), 'enable', 'off')
end
set(findobj(fig, 'tag', 'rating'),  'callback', ['web(''' selectedPlugins(lastSelect).webrating ''', ''-browser'');' ]);

tmpobj.selection = listboxVal;
set(fig, 'userdata', tmpobj);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function vl_simplenn_diagnose(net, res)
% VL_SIMPLENN_DIAGNOSE  Plot diagnostic information
%   VL_SIMPLENN_DIAGNOSE(NET, RES) plots in the current window
%   the average, maximum, and miminum element for all the filters
%   and biases in the network NET. If RES is also provided, it will
%   plot the average, minimum, and maximum element for all the
%   intermediate responses and deriviatives stored in RES as well.
%
%   This function can be used to rapidly glance at the evolution
%   of the paramters during training.

n = numel(net.layers) ;
fmu = NaN + zeros(1, n) ;
fmi = fmu ;
fmx = fmu ;
bmu = fmu ;
bmi = fmu ;
bmx = fmu ;
xmu = fmu ;
xmi = fmi ;
xmx = fmx ;
dxmu = fmu ;
dxmi = fmi ;
dxmx = fmx ;
dfmu = fmu ;
dfmi = fmu ;
dfmx = fmu ;
dbmu = fmu ;
dbmi = fmu ;
dbmx = fmu ;

for i=1:numel(net.layers)
  ly = net.layers{i} ;
  if ismember(ly.type, {'conv', 'bnorm'}) && numel(ly.filters) > 0
    x = gather(ly.filters) ;
    fmu(i) = mean(x(:)) ;
    fmi(i) = min(x(:)) ;
    fmx(i) = max(x(:)) ;
  end
  if ismember(ly.type, {'conv', 'bnorm'}) && numel(ly.biases) > 0
    x = gather(ly.biases) ;
    bmu(i) = mean(x(:)) ;
    bmi(i) = min(x(:)) ;
    bmx(i) = max(x(:)) ;
  end
  if nargin > 1
    if numel(res(i).x) > 1
      x = gather(res(i).x) ;
      xmu(i) = mean(x(:)) ;
      xmi(i) = min(x(:)) ;
      xmx(i) = max(x(:)) ;
    end
    if numel(res(i).dzdx) > 1
      x = gather(res(i).dzdx);
      dxmu(i) = mean(x(:)) ;
      dxmi(i) = min(x(:)) ;
      dxmx(i) = max(x(:)) ;
    end
    if ismember(ly.type, {'conv', 'bnorm'}) && numel(res(i).dzdw{1}) > 0
      x = gather(res(i).dzdw{1}) ;
      dfmu(i) = mean(x(:)) ;
      dfmi(i) = min(x(:)) ;
      dfmx(i) = max(x(:)) ;
    end
    if ismember(ly.type, {'conv', 'bnorm'}) && numel(res(i).dzdw{2}) > 0
      x = gather(res(i).dzdw{2}) ;
      dbmu(i) = mean(x(:)) ;
      dbmi(i) = min(x(:)) ;
      dbmx(i) = max(x(:)) ;
    end
  end
end

if nargin > 1
  np = 6 ;
else
  np = 2 ;
end

clf ; subplot(np,1,1) ;
errorbar(1:n, fmu, fmi, fmx, 'bo') ;
grid on ;
xlabel('layer') ;
ylabel('filters') ;
title('coefficient ranges') ;

subplot(np,1,2) ;
errorbar(1:n, bmu, bmi, bmx, 'bo') ;
grid on ;
xlabel('layer') ;
ylabel('biases') ;

if nargin > 1
  subplot(np,1,3) ;
  errorbar(1:n, xmu, xmi, xmx, 'bo') ;
  grid on ;
  xlabel('layer') ;
  ylabel('x') ;

  subplot(np,1,4) ;
  errorbar(1:n, dxmu, dxmi, dxmx, 'bo') ;
  grid on ;
  xlabel('layer') ;
  ylabel('dzdx') ;

  subplot(np,1,5) ;
  errorbar(1:n, dfmu, dfmi, dfmx, 'bo') ;
  grid on ;
  xlabel('layer') ;
  ylabel('dfilters') ;

  subplot(np,1,6) ;
  errorbar(1:n, dbmu, dbmi, dbmx, 'bo') ;
  grid on ;
  xlabel('layer') ;
  ylabel('dbiases') ;
end

drawnow ;

author: huangzehao

repo: caffe-vdsr


text:
function [pass, u1, u2, info1, info2] = test_scalarODE_sign(pref)
% A nonlinear CHEBOP test. This test tests a scalar ODE, where there is a
% breakpoint in the domain of the CHEBOP. Furthermore, the operator has a
% discontinuous coefficient, which will induce a further breakpoint in the
% solution. The problem is solved using chebcolloc1, chebcolloc2 and ultraS
% discretizations. The problem solved does not require damping for the Newton
% iteration to converge.
%
% Asgeir Birkisson, May 2014.

%% Setup
if ( nargin == 0 )
    pref = cheboppref;
end
dom = [-1 .5 1];

N = chebop(@(x,u) diff(u,2) + sign(x).*sin(u), dom);
N.lbc = @(u) u - 2;
N.rbc = @(u) u - 2;
rhs = 0;

%% Try different discretizations
% Start with chebcolloc2
pref.discretization = @chebcolloc2;
[u1, info1] = solvebvp(N, rhs, pref);

%% Change to ultraS
pref.discretization = @ultraS;
[u2, info2] = solvebvp(N, rhs, pref);

%% Change to chebcolloc1
pref.discretization = @chebcolloc1;
[u3, info3] = solvebvp(N, rhs, pref);

%% Did we pass?
% To pass, both residuals have to be small, but we should not expect u1 and u2
% to be identical!

% TODO: This used to be pref.bvpTol. Once we tune the algorithms better, should
% try to restore it.
tol = 1e3*pref.bvpTol;
pass(1) = norm(N(u1)) < tol;
pass(2) = norm(N(u2)) < tol;
pass(3) = norm(N(u3)) < tol;
pass(4) = ( (norm(u1 - u2) ~= 0) && (norm(u2 - u3) ~= 0) && ...
    (norm(u1 - u3) ~= 0));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function aux_SmrCumT(params, hParentFigure)
    % Plots seismic moment release and number of events versus time
    % function aux_SmrCumT(params, hParentFigure)
    %
    %
    % Incoming variables:
    % params        : all variables
    % hParentFigure : Handle of the parent figure
    %
    % J.Woessner, jowoe@gps.caltech.edu
    % updated: 09.11.05

    % Get the axes handle of the plotwindow
    axes(sv_result('GetAxesHandle', hParentFigure, [], guidata(hParentFigure)));
    set(gca,'NextPlot','add');
    % Select a point in the plot window with the mouse
    [fX, fY] = ginput(1);
    fprintf('X: %g Y: %g\n', fX, fY);
    % Plot a small circle at the chosen place
    plot(fX,fY,'ok');

    % Get closest gridnode for the chosen point on the map
    [fXGridNode fYGridNode,  nNodeGridPoint] = calc_ClosestGridNode(params.mPolygon, fX, fY);
    plot(fXGridNode, fYGridNode, '*r');
    set(gca,'NextPlot','replace');

    % Get the data for the grid node
    mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNodeGridPoint}, :);

    if (params.nGriddingMode == 1 && params.bMap)
        vDistances_ = sqrt(((mNodeCatalog_(:,1)-fXGridNode)*cosd(fYGridNode)*111).^2 + ((mNodeCatalog_(:,2)-fYGridNode)*111).^2);
        vSel = (vDistances_ <= params.fRadius);
        fCheckDist = max(vDistances_(vSel, :));
        mNodeCatalog_=mNodeCatalog_(vSel,:);
    elseif (params.nGriddingMode == 1 && ~params.bMap)
        [nRow_, nColumn_] = size(mNodeCatalog_);
        vXSecX_ = mNodeCatalog_(:,nColumn_);  % length along x-section
        vXSecY_ = (-1) * mNodeCatalog_(:,7);
        vDistances_ = sqrt(((vXSecX_ - fXGridNode)).^2 + ((vXSecY_ - fYGridNode)).^2);
        vSel = (vDistances_ <= params.fRadius);
        mNodeCatalog_=mNodeCatalog_(vSel,:);
        fCheckDist = max(vDistances_);
    end

    % Sort the catalog according to time
    [s,is] = sort(mNodeCatalog_(:,3));
    mNodeCatalog_ = mNodeCatalog_(is(:,1),:);
    % Calculate moment release
    [fCumMoment, vCumMoment, vMoment] = calc_moment(mNodeCatalog_);

    if exist('new_fig','var') &&  ishandle(new_fig)
        set(0,'Currentfigure',new_fig);
        disp('Figure exists');
    else
        new_fig=figure_w_normalized_uicontrolunits('tag','bnew','Name','Cumulative FMD and b-value fit','Units','normalized','Nextplot','add','Numbertitle','off');
        new_axs=axes('tag','ax_bnew','Nextplot','add','box','on');
    end
    set(gca,'tag','ax_bnew','Nextplot','replace','box','on');
    axs5=findobj('tag','ax_bnew');
    axes(axs5(1));
    vTime = mNodeCatalog_(:,3);
    vCumNumber = (1:length(vTime));
    vCumNumber(length(vTime)) = length(vTime);
    %plot(mNodeCatalog_(:,3),vCumMoment,'Linewidth',2,'Color',[0.7 0.7 0.7])

    [AX,H1,H2] = plotyy(mNodeCatalog_(:,3),vCumNumber,mNodeCatalog_(:,3),vCumMoment,'plot');
    % Figure refinements
    set(gca,'LineWidth',2,'FontSize',12','FontWeight','bold')
    set(AX(2),'YColor','k','XTicklabel',[])
    set(get(AX(1),'Ylabel'),'String','Number of events','FontWeight','bold','FontSize',12,'Color',[0 0 0]);
    set(get(AX(2),'Ylabel'),'String','Cumulative Moment [Nm]','FontWeight','bold','FontSize',12,'Color',[0.8 0 0]);
    set(H1,'Color',[0 0 0],'Linewidth',2);
    set(H2,'Color',[0.8 0 0],'Linewidth',2);

    xlabel('Time [years] ','FontWeight','bold','FontSize',12);
    set(AX(2),'XColor','k','YColor',[0.8 0 0])
    set(AX(1),'YColor','k')
end

author: CelsoReyes


text:
function [initState, transmat, mu, Sigma] = gausshmm_train_observed(obsData, hiddenData, ...
						  nstates, varargin)
% GAUSSHMM_TRAIN_OBSERVED  Estimate params of HMM with Gaussian output from fully observed sequences
% [initState, transmat, mu, Sigma] = gausshmm_train_observed(obsData, hiddenData, nstates,...)
%
% INPUT
% If all sequences have the same length
% obsData(:,t,ex)
% hiddenData(ex,t)  - must be ROW vector if only one sequence
% If sequences have different lengths, we use cell arrays
% obsData{ex}(:,t)
% hiddenData{ex}(t)
%
% Optional argumnets
% dirichletPriorWeight - for smoothing transition matrix counts
%
% Optional parameters from mixgauss_Mstep:
% 'cov_type' - 'full', 'diag' or 'spherical' ['full']
% 'tied_cov' - 1 (Sigma) or 0 (Sigma_i) [0]
% 'clamped_cov' - pass in clamped value, or [] if unclamped [ [] ]
% 'clamped_mean' - pass in clamped value, or [] if unclamped [ [] ]
% 'cov_prior' - Lambda_i, added to YY(:,:,i) [0.01*eye(d,d,Q)]
%
% Output
% mu(:,q)
% Sigma(:,:,q)

[dirichletPriorWeight, other] = process_options(...
    varargin, 'dirichletPriorWeight', 0);

[transmat, initState] = transmat_train_observed(hiddenData, nstates, ...
						'dirichletPriorWeight', dirichletPriorWeight);

% convert to obsData(:,t*nex)
if ~iscell(obsData)
  [D T Nex] = size(obsData);
  obsData = reshape(obsData, D, T*Nex);
else
  obsData = cat(2, obsData{:});
  hiddenData = cat(2,hiddenData{:});
end
[mu, Sigma] = condgaussTrainObserved(obsData, hiddenData(:), nstates, varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function imgRec = BoschettiDec(name)
%
%
%       imgRec = BoschettiDec(name)
%
%       Input:
%           -name: the prefix of the compressed HDR images using Boschetti
%           et al. method.
%
%
%     Copyright (C) 2012  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

%Read metadata
info = imfinfo([name, '_bos_RGB.jp2']);
decoded = sscanf(cell2mat(info.Comments), '%g', 3);
nBit = decoded(1);
maxE = decoded(2);
minE = decoded(3);

maxVal = 2^nBit - 1;

%Reading and Decoding Eq
EqDec = double(imread([name, '_bos_E.jp2'])) / maxVal;
EDec = EqDec * (maxE - minE) + minE;
mult = 2.^EDec;

%Decoding RGB
RGBDec = double(imread([name, '_bos_RGB.jp2'])) / maxVal;

%Reconstruction
imgRec = zeros(size(RGBDec));
for i=1:size(imgRec, 3)
    imgRec(:,:,i) = (RGBDec(:,:,i)) .* mult;
end

end

author: banterle


text:
function dif = r8mat_diff_frobenius ( m, n, a, b )

%*****************************************************************************80
%
%% R8MAT_DIFF_FROBENIUS: Frobenius norm of the difference of two R8MAT's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, real A(M,N), B(M,N), the matrices for which we
%    are to compute the Frobenius norm of the difference.
%
%    Output, real DIF, the Frobenius norm of A-B.
%
  dif = sqrt ( sum ( sum ( ( a(1:m,1:n) - b(1:m,1:n) ).^2 ) ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_failed = test_libltfat_fft(varargin)
test_failed = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

[~,~,enuminfo]=libltfatprotofile;
phaseconv = enuminfo.ltfat_phaseconvention;

fftwflags = struct('FFTW_MEASURE',0,'FFTW_ESTIMATE',64,'FFTW_PATIENT',32,'FFTW_DESTROY_INPUT',1,...
    'FFTW_UNALIGNED',2,'FFTW_EXHAUSTIVE',8,'FFTW_PRESERVE_INPUT',16);

Larr  = [351 350   9   1];
Warr  = [  10       3   3   1];

for idx = 1:numel(Larr)
    L = Larr(idx);
    W = Warr(idx);

    f = randn(L,W,flags.complexity) + 1i*randn(L,W,flags.complexity);
    fin = complex2interleaved(f);
    fPtr = libpointer(dataPtr,fin);

    c = cast(randn(L,W)+1i*randn(L,W),flags.complexity);
    cout = complex2interleaved(c);
    coutPtr = libpointer(dataPtr,cout);

    truec = fft(f);

    funname = makelibraryname('fft',flags.complexity,0);
    status = calllib('libltfat',funname,fPtr,L,W,coutPtr);

    res = norm(truec - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
    fprintf(['FFT   L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

    % With plan
    c = cast(randn(L,W)+1i*randn(L,W),flags.complexity);
    cout = complex2interleaved(c);
    coutPtr = libpointer(dataPtr,cout);

    plan = libpointer();
    funname = makelibraryname('fft_init',flags.complexity,0);
    statusInit = calllib('libltfat',funname,L,W,fPtr,coutPtr, fftwflags.FFTW_MEASURE, plan);

    funname = makelibraryname('fft_execute',flags.complexity,0);

    statusExecute = calllib('libltfat',funname,plan);

    res = norm(truec - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['FFT L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

    c = cast(randn(L,W)+1i*randn(L,W),flags.complexity);
    cout = complex2interleaved(c);
    coutPtr = libpointer(dataPtr,cout);

    funname = makelibraryname('fft_execute_newarray',flags.complexity,0);
    statusExecute = calllib('libltfat',funname,plan,fPtr,coutPtr);

    res = norm(truec - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['FFT L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

    funname = makelibraryname('fft_done',flags.complexity,0);
    statusDone = calllib('libltfat',funname,plan);

    %%%%%% Inplace
    c = f;
    cout = complex2interleaved(c);
    coutPtr = libpointer(dataPtr,cout);

    plan = libpointer();
    funname = makelibraryname('fft_init',flags.complexity,0);
    statusInit = calllib('libltfat',funname,L,W, coutPtr, coutPtr, fftwflags.FFTW_MEASURE, plan);

    funname = makelibraryname('fft_execute',flags.complexity,0);
    statusExecute = calllib('libltfat',funname,plan);

    res = norm(truec - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['FFT L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

    c = f;
    cout = complex2interleaved(c);
    coutPtr = libpointer(dataPtr,cout);

    funname = makelibraryname('fft_execute_newarray',flags.complexity,0);
    statusExecute = calllib('libltfat',funname,plan,coutPtr,coutPtr);

    res = norm(truec - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['FFT L:%3i, W:%3i, %s %s %s\n'],L,W,flags.complexity,ltfatstatusstring(status),fail);

    funname = makelibraryname('fft_done',flags.complexity,0);
    statusDone = calllib('libltfat',funname,plan);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function oR = cleanUp(oR)

% testing:
% cs = crystalSymmetry('222');
% oR = cs.fundamentalRegion;

% we need to restrict to the relevant normals
% this works nice if rotations about 180 degree are present
is180 = abs(oR.N.angle)>pi - 1e-3;
if ~any(is180)

  % then the normals should be inside itself
  Nq = unique(oR.N);
  Nq = Nq.subSet(oR.checkInside(-reciprocal(Nq)));
  oR.N = Nq;

else % otherwise the condition seems to strict

  % we consider first the case without 180 degree rotations
  oRfull = oR;
  oRfull.N(is180) = [];
  oRfull = oRfull.cleanUp;

  ind = oR.checkInside(oRfull.V);

  inside = cellfun(@(x) any(ind(x)), oRfull.F);

  Nq = [oRfull.N(inside),oR.N(is180)];
  oR.N = Nq;
end

nqa = unique(Nq.axis,'antipodal');

if isempty(nqa)

  oR.V = orientation(oR.CS1,oR.CS2);
  oR.F = {};
  return

elseif length(nqa)==1

  v = perp(nqa);
  rot = rotation.byAxisAngle(nqa,(0:90:270)*degree);
  oR.V = orientation.byAxisAngle(rot*v,pi*ones(4,1),oR.CS1,oR.CS2);
  oR.F = repcell((1:4).',size(oR.N));

  return
end

% compute vertices

% set up faces
for j = 1:length(Nq)

  Nqj = Nq(j);
  oNq = Nq; oNq(j) = [];
  oNq(abs(dot(oNq.axis,Nqj.axis))>1-1e-6) = [];

  % find closest other normal
  [~,j0] = min(angle(Nqj,oNq) + angle(Nqj.axis,oNq.axis));

  % rotate axis such that current face is z-direction
  % and closest plane
  aNq = oNq.axis;
  aNq0 = aNq(j0);
  aNqj = Nqj.axis;
  aNq0 = aNq0 - dot(aNq0,aNqj) * aNqj;
  rot = rotation.map(aNqj,zvector,aNq0,xvector);
  aNq = rot * aNq; % the rotated axes

  % order the other normals according to
  % 1. angle to aNq0
  % 2. angle to Nqj
  [~,order] = sort(100*mod(aNq.rho+1e-5,2*pi) + aNq.theta);
  order(end+1) = order(1); %#ok<AGROW>

  % the first edge
  e = order(1);
  % go through all potential edges and compute vertices
  V{j} = quaternion;
  for i = 2:length(order)

    % if axis is in the same direction -> skip
    if isnull(aNq(order(i)).rho - aNq(e).rho), continue; end

    % is axis is inline whith Nqj -> skip
    %if abs(dot(aNq(order(i)),zvector))>1-1e-6, continue; end

    % compute vertice
    v = cross(Nqj,oNq(e),oNq(order(i)));
    v = v ./ norm(v);
    if abs(angle(v * Nqj))<1e-5 , continue; end

    % check vertex is inside
    if ~oR.checkInside(v), continue; end

    % if the rotational angle is 180 degree then +- axis is possible
    % take the quaternion is the axis inside the fundamental region

    if abs(v.a)<1e-5, v.a = sign(v.a)*1e-5; end

    if ~oR.axisSector.checkInside(v.axis) && (abs(v.a) < 1e-4)
      v.a = -v.a;
    end

    % if we have found a vertice - store it
    V{j}(end+1) = v./norm(v);
    e = order(i);

  end
end

% make a unique list of vertices
sV = cellfun(@length,V);
[V,~,iV] = unique([V{:}]);
oR.V = orientation(V,oR.CS1,oR.CS2);
sV = sV(sV>0);
F = 1:sum(sV);
F = mat2cell(F,1,sV);
F = cellfun(@(x) iV(x),F,'uniformOutput',false);
oR.F = cellfun(@(x) x([diff(x)~=0; true]),F,'uniformOutput',false);

end

function q = reciprocal(q)

q = q .* rotation.byAxisAngle(-q.axis,pi);

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
%  INTERNAL FUNCTION: Projection of covariance matrix such the eigenvalues
%  are sandwiched.
%
%  ::
%
%     vcov = project(vcov0,e_min,e_max);
%
%  Args:
%
%     - **vcov0** [matrix]: initial covariance matrix
%     - **e_min** [[]|{sqrt(eps)}]: scalar such that the minimum eigenvalue of vcov
%       is greater than or equal to "e_min".
%     - **e_max** [[]|{1/e_min}]: scalar such that maximum eigenvalue of vcov
%       is less than or equal to "e_max"
%
%  Returns:
%     :
%
%     - **vcov** [matrix]: updated covariance matrix
%
%

author: jmaih